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ABSTRACT 

We present a semi-analytic model for the interstellar medium that considers local 
processes and structures of turbulent star-forming gas. A volume element of the inter- 
stellar medium is described as a multi-phase system, comprising a cold and a warm 
gas phase in effective (thermal plus turbulent) pressure equilibrium, and a stellar 
component. The cooling instability of the warm gas feeds the cold phase, while var- 
ious heating processes transfer cold gas to the warm phase. The cold phase consists 
of clumps embedded in diffuse warm gas, where only the molecular fraction of the 
cold gas may be converted into stars. The fraction of molecular gas is approximately 
calculated, using a Stromgren-like approach, and the efficiency of star formation is 
determined by the state of the cold gas and by the turbulent velocity dispersion on 
the clump length scale. Gas can be heated by supernovae and UV-emission of massive 
stars, according to the evolutionary stages of the stellar populations and the initial 
mass function. Since turbulence has a critical impact on the shape of the gaseous 
phases, on the production of molecular hydrogen and on the formation of stars, the 
consistent treatment of turbulent energy - the kinetic energy of unresolved motions 

- is an important new feature of our model. Besides turbulence production by super- 
novae and by the cooling instability, we also take into account the forcing by large 
scale motions. 

We formulate a set of ordinary differential equations, which statistically describes 
star formation and the exchange between the different budgets of mass and energy 
in a region of the interstellar medium with given mean density, size, metallicity and 
external turbulence forcing. By exploring the behaviour of the solutions, we find equi- 
librium states, in which the star formation efficiencies are consistent with observations. 
Kennicutt-Schmidt-like relations naturally arise from the equilibrium solutions, while 
conventional star formation models in numerical simulations impose such relations 
with observed efficiency parameters as phenomenological calibrations. 

Beyond the semi-analytic approach, a potential application is a complete subgrid 
scale model of the unresolved multi-phase structure, star formation and turbulence in 
simulations of galaxies or in cosmological simulations. The formulation presented in 
this article combines various models focusing on particular processes and yet can be 
adopted to specific applications, depending on the range of resolved length scales. 
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- turbulence 



1 INTRODUCTION 

The capabilities of contemporary supercomputing enable 
us to model the evolution of the baryonic gas in the uni- 
verse with unprecedented sophistication. Adaptive meth- 
ods such as smoothed particle hydrodynamics (SPH) and 
adaptive mesh refinement (AMR) in Eulerian grid codes al- 
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low us to cover a huge dynamic range such that simula- 
tions of the formation and evolution of galaxies from cosmo- 
logical initial conditions at high resolution (~ 100 pc) ar e 
within reach |Gnedin fc KravtsovllioTol : I Agertz et al]|201ll ). 
In simulations of isolated disc galaxies, it is feasible to re- 
solve length scales down to ~ 10 pc |Agertz et al.l 120091 : 
iTasker fc Tanll2009h . Computations on these length scales 
entail the problem to account for variou s physical processe s 
in the multi-phase interstellar medium (|Maver et al. 2008). 
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Notwithstanding the high numerical resolution that can be 
achieved, several important processes cannot be fully re- 
solved and have to be described by means of a sub-grid scale 
(SGS) model. 

The distribution of the gas among the different phases 
of the ISM is controlled by the following physical pro- 
cesses. The fragmentation of warm neutral gas (num- 
ber density n < 1cm -3 , temperature T > 10 4 K) is 
driven by gravi tational instabilities on length scales ~ 
0.1 lkpc (e.g. lToomrdll964l;IWada et al.ll2002l; iKravtsov 
20031 : iLi et al.ll2005l : IWada fc Normanl 120071 : lAgertz et al. 
20091 ). The gravitational contraction of gas is supported 
by cooling processes in converging flows that produce the 
cold neutral phase (n > 10cm~ 3 , T < 10 3 K). Gravi- 
tational and cooling instabilities, and, possibly, magnetic 
fields act in concert to form dense star-forming clouds, 
in which molec ular hydrogen is p r oduced at densities > 
100 cm" 3 (e.g. lDobbset al.1 120081 : iRobertson fc Kravtsovl 
120081 ; iTasker fc Tarj|200d ). The radiation from hot massive 
stars and blast waves from supernovae feeds energy back into 
the interstellar medium. Gravity, cooli ng, and stellar feed- 
back are potential drivers of turbulence jElmegreen fc Scalo 
20041 : iMac Low fc Klessenll2004l:lde Avillez fc Breitschwerdt 
2004: iBurkert et alJ l20ld: iKlessen fc Hennebelld l20ld ; 



proposed (e . g., IGnedin et alj 20091: Krumholz et al. 2009 



zuum; icurKert et ai. izuiu: rvicsscn ncnncp ena izuiu 
Bournaud et al.ll20ld : lFederrath et al which, in turn 



has an impact on the st ability of the gas i Bonazzola et all 
1 19921 : iRomeo et al1l2010j ) . 

In large-scale simulations, where the smallest resolved 
length scales range from the scale of star-forming regions 
to galactic scales, it is a major challenge to account for 
the sub-resolution structure and dynamics of the ISM (re- 
cent reviews are given by iMcKee fc Ostriken (20071 ) and 
iHenslerl (200d )). On the one hand, isolated disc galaxy sim- 
ulations serve as idealized models of galaxy evolution that 
avoid so me of the difficulties one faces in cosmological sim- 
ulations JPobbs et all 120081: IRobertson fc Kravtsovl |200S|; 



Agertz et all 120091: ITasker fc Tanl 120091 : iDobbs fc Pringld 



20ld : iBournaud et al.ll20ld ). Although the highly artificial 



initial conditions are problematic, isolated discs can be used 
to study dynamical properties of the ISM at high resolu- 
tions and to test advanced models of unresolved processes. 
Because artefa cts may result from discs that are adiabati- 
cally unstable, IWang et all (|2010h defined an adiabatic disc 
that is stable over the rotation time scale. On the other 
hand, substantial efforts have been made to zoom into ha- 
los from cosmological simulations and to re-simulate galax- 
ies evol ving from those halos a t the highest feasib l e res- 
olution jGoyernato et al. 2007 1_ Gnedin fc Kravtsovl 2010l : 



lAgertz et al.ll201ll : lGovernato et al.ll2010l : iGreif et aljkoid ) 

Several models were developed i n the past to de 



scibe the multi-phase I S M (e. 



19771: lYepes et all Il997l: iGnedhf Il998l: iKlvpinl 1 19981 : 
Hultman fc Pharasvnl 1 19991 : IStinson et al.ll200d ). An often 
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used type of model for star formation and stellar feedback in 
cosmological smoothe d particle hydrodynam i cs (SP H) simu- 
lations is described in lSpringel fc Hernquistl d2003l) fSH03), 
which is an adaption of the model introduced bv lYepes et alJ 
l| 19971 ). Basically, rate equations for the densities of the cold 
and hot gas phases are formulated, including sources and 
sinks related to star formation and feedback from super- 
novae. Recently, a variety of phenomenological models that 
treat particular physical processes in the ISM have been 



[ Joung et aill2009l:lMurante et alj|20ld ; lOstriker et allboiol : 
IPadoan fc Nordlundll201ll ). Some of these models are de- 
signed to account for sub-grid scale physics in numerical 
simulations. Others are mainly intended to obtain analyti- 
cal or semi-analytical predictions that can be compared to 
observations. Even so, particular components of the latter 
class of models could be incorporated into an SGS model. 
In the following, we briefly review these models from the 
perspective of the physical processes involved. 

IPadoan fc Nordlundl (201ll ) [PN11] parametrize the star 
formation rate per free-fall time as a function of the virial 
parameter, i. e., the turbulent velocity dispersion relative 
to the specific gravitational energy, by using data from 
forced isothermal MHD tu rbulence simulations. Following 
iKrumholz fc McKed (20051 ) [KM05], the star formation rate 
is calculated by integrating density fluctuations beyond a 
critical density that is given by the virial parameter and the 
Mach number of the tur bulent cold neutral med ium. How- 
ever, as pointed out by IKrumholz et all (2009) [KMT09], 
new observations reveal a tight correlation between the 
molecular hydrogen surface density and the star formation 
rate. They present an analytic model that includes approx- 
imate calculations of molecular hydrogen fraction from a 
spherical-cloud model and the star formation efficiency per 
free-fall time on the basis of the numerical parametrization 
in KM05. This model reproduces the Kennicutt-Schmidt re- 
lation between the star formation rate and the surface den- 
sity on length scales of the order of a kpc in recent surveys. 

By assuming a constant star formation efficiency, the 
formation of molecular hydrogen in cosmological simula- 
tions is modelled by an ap proximate treatment of shield- 
ing and photo-dissociation in IGnedin et all (2009) [GTK09]. 
As in KMT09, the star formation rate is assumed to 
be proportional to the molecular hydrogen density rather 
than the density of the cold neutral medium. The un- 
resolved density structure of the gas is parametrized by 
a clumping factor, and the efficiency of star formation 
per free-fall ti me in molecular clouds is set to 1 %. Us- 
ing this model, IGnedin fc Kravtsovl (20ld ) investigate the 
Kennicutt-Schmidt relation in galaxies at high redshifts. For 
simulati ons of isolated discs with molecular hydrogen chem - 
istry, see lDobbs et al. I (20081) : IRobertson fc Kravtsovl (20081 ). 

The KMT09 and GTK09 models focus on molecular hy- 
drogen to predict the star formation rate, whereas the multi- 
phase structure and the turbulent dy namics of the I S M are 
not addressed explicitly. In contrast, iKoppen et all (19981 ) 
formulate a dynamical model for the evolution of a massive 
and a low-mass star component and clouds embedded in hot 
gas, with various interaction pr o cesses . In a similar way, the 
model of lSpringel fc Hernquistl (2003h considers interacting 
cold and warm phases and stars. A simple multi-phase SGS 
m odel of star formation and supernova feedback is proposed 
bv lMurante et all (20ld ). By assuming that the amount of 
molecular hydrogen is controlled by the pressure of the ISM, 
rate equations for the mass and the energy of a cold and 
a warm phase are so lved in addition to t he mass that is 
converted into stars. lOstriker et"aH (20ld ) present a con- 
siderably more detailed analytical model that separates the 
ISM into a diffusive gas component and into gravitationally- 
bound clouds, in which stars are formed at a rate that is pro- 
portional to their mass. The basic parameters of this model 
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are the ratio of the thermal to the effective pressure (i. e., 
the sum of thermal, turbulent and magnetic pressures) of 
the diffusive gas, the fraction of warm diffusive gas (com- 
plementing a cold diffusive phase), and the star formation 
efficiency of the clouds. The main idea is that the radiation 
of young massive stars heats the diffusive ISM and the mass 
exchange between the diffusive components and the clouds 
regulates star formation. Turbulence and the conversion of 
atomic into molecular hydrogen are not decisive for the reg- 
ulation process. 

To i nclude SN feedback i n cosmological simulations, for 
example. IStinson et al.l l|200rj ) model the impact of SN blast 
waves on the ther mal structure of the ISM. In contrast, 
iJoung et all ((2009) [JMB09] propose a non-thermal treat- 
ment of SN feedback. They formulate a dynamical equation 
to compute the numerically unresolved turbulent pressure of 
the ISM, with the rate of energy injection by SN blast waves 
as source term (internal turbulence driving). The turbulent 
pressure is proportional to the energy density of numerically 
unresolved turbulent velocity fluctuations. The coefficients 
of the equation for the turbulent pressure are calibrated on 
galactic-sc ale simulations of the ISM. A s imilar approach is 
utilised in IScannapieco fc Briiggenl (|2010| ) [SB10] for turbu- 
lence in galaxy outflows. 

Although feedback models using the turbulent pressure 
are promising, JMB09 and SB 10 do not account for the in- 
crease of the turbulent pressure by the energy transfer from 
resolved to unresolved scales via the turbulent cascade (ex- 
ternal turbulence driving). We expect this production chan- 
nel to be important because turbulence in the ISM is to 
some extent dri ven by gravitationa l instabilities on large, 
galac tic scales (jRomeo et alJ l2010l : iKlessen fc Hennebelld 
20ld). For the local comput ation in numerical simulations, 
Schmidt fc Federratbl (j201lh [SF11] formulated and tested 
an SGS model for highly compressible turbulence. This 
model is also based on a dynamical equation for the nu- 
merically unresolved turbulent energy. However, in addition 
to diffusion and dissipation terms, SGS turbulent energy is 
produced by the shear of resolved small-scale fluctuations, 
i. e., the turbulent cascade. The rate of production by the 
turbulent cascade is called the turbule nt energy flux. Simu- 
lations of forced super sonic turbulence ( Schmidt et al.ll2009l ; 
iFederrath et al.ll2010h were used to verify a new closure for 
the compressible turbulent energy flux. In large eddy simu- 
lations (LES), a closure is an approximation to a quantity 
that depends on unresolved density and velocity fluctua- 
tions. Moreover, it is demonstrated that the SGS model ful- 
fils several basic requirements, such as a constant mean dis- 
sipation rate, independent of the numerical resolution, and 
a power-law scaling of the SGS turbulent energy. For com- 
pressible turbulence driven by large-scale instabilities, this 
SGS model is the only model for computing the turbulent 
pressure consistently that has been systematically tested so 
far. Feedback can be included as an additional production 
term in the SGS turbulent energy equation. Since the unre- 
solved turbulent velocity fluctuations in galaxy simulations 
are comparable to the speed of sound, we expect significant 
effects of the corresponding turbulent pressure, particularly 
with regard to the regulation of star formation. 

The aim of this work is to bring together different ap- 
proaches, using the SH03 model as a basic framework. Our 
treatment of star formation and molecular hydrogen for- 



mation is guided by KM05, KMT09, and PN11. To heat 
the interstellar gas, Lyman-continuum radiation of young 
massive stars and supernova feedback are calcul ated from 
the m odelled star formation history, assuming the lChabrierl 
(2001) initial mass function. We incorporate internal tur- 
bulence driving by the thermal instability and by a non- 
thermal fraction of supernovae feedback, as in JMB09. By 
adding the turbulent pressure to the thermal pressure, tur- 
bulence influences the pressure balance between the phases 
and, in the highly turbulent regime, it significantly affects 
the gravitationally unstable mass fraction in the cold-gas 
phase. The key to the fluid-dynamical computation of the 
turbulent energy, including external driving via a turbulent 
cascade, is SF11. 

In this paper, we devise a semi-analytic formulation to 
describe the evolution of the two gas phases, turbulence, 
star formation, and feedback by averaged quantities in a 
box of given size. These one-zone calculations allow us to 
investigate the dependence on the control parameters (total 
gas density, metallicity, constant rate of turbulent energy 
production by external driving) and the coefficients of the 
models. In particular, we calculate the star formation effi- 
ciency for self-regulated equilibria. These equilibrium solu- 
tions are useful in their own right for a parametrization of 
the star formation efficiency in various astrophysical appli- 
cations. The full implementation as a sub-grid scale model 
for cosmological and galaxy-scale simulations is the goal of 
future work. 

An outline of the proposed multi-phase model will be 
given in Sect. [5] followed by detailed descriptions of the 
star formation model (Sect. [3} and the model equations 
for the mass and energy budgets of the warm and cold 
phases (Sect. 3}. In Sect.0 we consider limiting cases (single 
phase, constant star formation rate in equilibrium). To test 
our model, we discuss results from one-zone calculations in 
Sect. [U including a comparison with observations. Finally, 
we present our conclusions and an outlook to the application 
of the model in numerical simulations. 



2 OUTLINE OF THE MODEL 

The base concept of this model is to split the density in a 
reference volume V — I 3 (i.e. a grid cell) into a cold and a 
warm pha se density with separa t e the rmal energy budgets, 
as used bv lSpringel fc Hernqui st (2003). The separation into 
two phases results from the cooling instability. In addition 
to the thermal energies of the cold and warm gas, the tur- 
bulent energy on the length scale I is computed. Contrary to 
most star formation models that are used in contemporary 
numerical simulations, we determine the star formation effi- 
ciency per free-fall time scale based on local properties and 
processes of the turbulent multi-phase medium. To calculate 
the star formation efficiency, the typical length scale of cold- 
gas clumps embedded in the warm neutral medium and the 
fraction of molecular hydrogen are important parameters. 
The molecular hydrogen fraction, in turn, depends on the 
composition and the density of the gas. To close the system 
of equations, we assume virial equilibrium for the cold phase, 
which is largely dependent on the effective pressures, i. e., 
the sum of thermal and turbulent pressures, of the phases. 



4 H. Brawn & W. Schmidt 

Table 1. Set of model parameters and important variables. 



Symbol 


Description 


main parameters 


I 


size of region 


P 


total mass density 


u c 


specific thermal energy of cold gas 


S 


rate of energy injection by the turbulent cascade 


l» 


intensity of incident UV-radiation field 


process parameters 




efficiency of cold phase evaporation by clump collisions 


ctt 


efficiency of turbulence production via phase separation 


^SN 


efficiency of turbulence production by SNe 


"SN 


specific energy of SN-ejecta 


V 


turbulent velocity scaling exponent 


b 


compressive factor, describing the ratio of solenoidal and compressive turbulent modes 


/loss 


fraction of mass ejected during prestellar collapse 


Cm 


fraction of newly build up metals in SN-ejecta 


£Lyc 


energy deposited in gas per absorbed Lyman continuum photon 


important variables 


u w 


specific thermal energy of warm gas 


et 


specific turbulent energy 


Pw 


fractional density of warm gas 


pw.pa 


average density in the warm phase 


Pc 


fractional density of cold gas 


pc.pa 


average density in the cold phase 


Ps 


averaged stellar mass density 


/c,H 2 


mass fraction of shielded molecular gas in the cold phase 


Z C 


size of cold clumps 


SFR c ,ff 


faction of shielded molecular gas converted into stars per respective free fall time 




faction of total density converted into stars per respective free fall time 


z 


mass fraction of heavy elements 



Since the turbulent pressure contribution is scale-dependent, 
the equilibrium also depends on the clump length scale. 

In the following, quantities with subscript 'c' belong to 
the cold phase, those with V to the warm phase, those 
with V to the star formation and those without the latter 
subscripts denote quantities of all the gas in the reference 
volume. An overview of used model specific parameters and 
variables is given in table [T] 

2.1 Specific energy variables 

The total thermal energy density up can be expressed as 
sum of the thermal energies of the cold and warm phases: 

up = u c p c + it w p w , (1) 

where fractional densities p c and p w are given by the 
masses m w and m c in the warm and cold phases, respec- 
tively, divided by the reference volume V, and p = p c + p w 
is the total gas density. 

The specific thermal energy of the warm phase, u w is 
changed by radiative cooling and heating, the mixing of hot 
SN-ejecta and cold gas, and turbulent dissipative heating. 
On the other hand, we assume that u c , the specific ther- 
mal energy of the cold phase, has a constant value, corre- 
sponding to an average temperature T c — 50 K of the cold 
phase. Numerical simulations suggest that the isothermal 



approximation is reasonable for the cold phase of the in- 
terstellar medium, because most of the gas in the cold gas 
is situated close to the asymptotically isothermal branch of 
the equilibrium curve between radiative coolin g and heating 
jSeifried et al.ll201ll ; lAudit fc Hennebelldl2oloh . To preserve 
energy conservation in our model, we account for any heat- 
ing process that affects the cold gas by a transfer of a certain 
amount of cold gas to the warm phase. 

Apart from the thermal energy, we assume that the gas 
in both phases has a certain specific turbulent energy et that 
corresponds to nearly isotropic random motions on length 
scales smaller than the size I of the reference volume. An 
exact definition of e t will be given on the basis of a decom- 
position of the fluid-dynamical equation in scale space. 

2.2 Density variables and effective pressure of the 
gas phases 

Since each phase fills only a fraction of the total volume V , 
we define the average densities within the phases, p c , P a and 
Pw.paEl by the identities 

m c = Pc.paK = PcV, (2) 
m w = p w ,paVw= pwV. (3) 

1 Subscript 'pa' means 'phase average' 
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where V c is the volume occupied by the cold gas phase, 
arid V w = V-Vc- 

Quantities such as the star formation rate, the molecu- 
lar fraction in the cold phase and the cooling rate depend on 
Pw, P a and Pc.pa- To determine V c , it would be necessary to 
know the structure of the two-phase medium on length scales 
smaller than I. In principle, one could parametrize the cold 
gas fraction V c /V fr om small-scale simu lations of thermally 
bistable turbulence (Sei fried et all 1201 If ). However, because 
of the high sensitivity of the thermal instability on the en- 
vironment (boundary conditions, gas density, etc.), it is not 
obvious how to relate the parameters of such idealised simu- 
lations to the the local properties of a grid cell in large-scale 
simulations. 

A much simpler approach is to assume that the cold 
gas is nearly in viral equilibrium if turbulence is accounted 
for and that clouds of cold gas with a characteristic scale l c 
are embedded in the warm phase. The effect of turbulence 
can be described by an effective pressure that includes both 
microscopic (thermal) and macroscopic (non-thermal) mo- 
tions (a precise definition will be given below). For a spher- 
ical cloud of density p c , P a and size l c , the generalized virial 
theorem implies the equilibrium condition 



3P C . 



-Gp 



r 



3P W 



o. 



(4) 



where the effective pressure of the war m phase is subs tracted 
as external pressure (see Sect. 14.1 in lLeaueux|[2005h . Since 
the turbulent pressure depends on the length scale, P c ,cff 
and P w , c ff are also functions of l c . In principle, this equation 
could be used to determine the length scale l c . It turns out, 
however, that the resulting system of equations is generally 
not well posed, meaning that no solutions exist for regions in 
the parameter space that definitely could be swept through 
in numerical simulations. As a consequence, either the rela- 
tively simple model with a single, characteristic length scale 
l c has to be abandoned or the assumption of virial equilib- 
rium as formulated above has to be loosened. A multi-scale 
model might eventually result from recent theoretical de- 
velopments (P. Hennebelle, private communication). In this 
article, we choose the second option and investigate its con- 
sequences. Typically, structures satisfying Eq. Q are not 
gravitationally bound. The dominant contributions come 
from the effective pressure, and these structures are held 
together by the pressure that is exerted by the surround- 
ing warm gas. For this reason, the gravitational energy term 
can be neglected, and we obtain an approximate effective 
pressure balance: 



(•») 



On the average, the turbulent pressure significantly con- 
tributes to the support of the cold gas against gravity. In 
order to connect the properties of the cold phase to the star 
formation rate, we assume that localized regions exist in the 
cold phase, where weak turbulent pressure support persists 
over sufficiently long periods of time so that the gas can col- 
lapse. The existence of such regions is a consequence of the 
intermittency of turbulence. The critical size of these regions 
is roughly given by the thermal Jeans length, 



Aj.c = c< 



(V- 

V7W>c, P 8 



1/2 



7r(7 — 1)m c 

Gp c ,pa 



1/2 



(6) 



where c c = [7(7 — l)^] 1 / 2 is the speed of sound in the cold 
gas, 7 the polytropic equation of state parameter and G the 
gravitational constant. Thus, we define the length scale l c 
by 



Zc = A, 



(7) 



The effective pressure equilibrium ((SJ) and the length 
scale (J7)l really have a complementary meaning. While the 
former statistically accounts for the overall effect of turbu- 
lence, the latter specifies a typical size of locally collapsing 
structures in the cold phase. In a certain sense, this cor- 
responds to the fact that molecular clouds do not collec- 
tively collapse although their mass is much greater than the 
thermal Jeans mass, while gravitationally unstable cores are 
formed locally ()Mac Low fc KlessenlfeoO^ . 

From the effective pressure balance ([5]) between the 
phases follows the ratio 



Pc,pa 
Pw.pa 



0"c,eff 



(8) 



where <r WiC ff and <r CiC ff are functions of the internal en- 
ergies u w and u c , and the turbulent energy et. Combining 
Eq. ([2}[8]), we can express the phase densities and volumes in 
terms of the fractional densities and the specific pressures: 

PcV 



pc.pa — Tw Pw -j- Pc, Vc 



— r vj Pc + P« 



Vw = 



pwV 



(9) 
(10) 



Pw ~\~ T w p c 

Furthermore, we define the stellar density p s to be the 
the stellar mass within the reference volume V divided by 
that volume: 



Ps = 



V 



(11) 



Numerical simulations of forced turbulence in thermally 
bistable gas indicate that the specific turbulent energy is 
nearly isotrop i c and uniformly distributed among the phases 
ISeifried et al.l l|201ll ). Thus, the turbulent velocity dispersion 
within the cold gas can be related to the turbulent energy 
on the length scale I via the power law 

2*7 



(12) 



The scaling exponent r\ is constrained by 1/3 ^ 77 ^ 1/2, 
where the lower and upper bounds correspond to Kol- 
mogorov and Burgers scaling, respectively. This scaling law 
is consist ent with the observed <r c -scaling relation (see, for 
example, lLarson|[l98ll ). 

With the above definition, the effective pressure of the 
cold gas on the length scale Z c is given by (see SF11) 



Cc . 2 

h cr c 

7 



fc.eff = Pc,paO" c ,cff = Pc 

= (7-l)Pc,paMc(l+|M^) 



(13) 



where M c = \f?>a c jc c is the root mean square Mach num- 
ber of turbulence in the cold phase. The turbulent pressure 
Pw,e« is given by an analogous expression, with cold-phase 
quantities replaced by the corresponding quantities in the 
warm phase. With these definitions, the variables l c , p c , P a, 
and p w ,pa can be determined solving Eqs. ©, (JTOj and (0 
iteratively. 
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2.3 Gas composition variables 

The chemical composition of the gas is for simplicity, as 
we do not track individual species, identified by its mass 
fraction of heavy elements Z. Given the mass fraction of 
helium Yq at solar metallicity Zq and its primordial value 
Yprim) Y at metallicity Z is assumed to be Y — Y pl im + 
(Yq — Yp^m) Z. Then the total mass fraction of hydrogen X 
is given by X — 1 — Y — Z. If the the gas is neutral, but not 
molecular, which is approximately true in the warm phase 
and in the cold phase gas outside of molecular cores, the 
mean molecular weight fj, is given by 

(be/i) _1 = Xrrijj 1 + Yrriftl + Zm z ', (14) 

where mn and mne are the atomic masses of hydrogen and 
helium, respectively, and mz is the average atomic mass of 
the heavier elements. Within the molecular cores of the cold 
phase we assume the gas to be fully molecular. 



3 STAR FORMATION 

Following KMT09, cold gas is converted into stars at a rate 
that depends on the mass of molecular hydrogen in the refer- 
ence volume (mHj = / c ,h 2 PcV , where / c ,h 2 i s the molecular 
hydrogen fraction in the cold gas phase): 

(1 — /loss)/c,H 2 Pc 



Ps = 



ts 



We define the star formation time scale t s by 

t c ,a 



SFR c , ff ' 



(15) 



(16) 



where the free-fall time scale in the cold gas is given by the 
phase-average (not the fractional) density: 

3tt 



i 2 



32Gp c , P a 

and SFRc,fi is the dimensionless star formation rate per 
free fall time t c ,s- Not all the mass in collapsing prestel- 
lar cores eventually ends up in a stars. A fraction /i oss ~ 
0.5 . . . 0.7 of mass is ejected during prestellar collapse 
fe.g.lMatzner fc McKell2000l : iHennebelle fc Chabrieij|20Qgl ; 
IChabrier fc HennebellelfeOlOl ). We account for the mass ejec- 
tion by correcting the star formation rate by the factor 
(1 - / loss ) in Eq. J25J). 

To calculate SFR c> fr, KMT09 derive a parametrization 
in terms of the gas column density, which reproduces impor- 
tant observational results from recent high-resolution sur- 
veys. These data also imply that the star formation is tightly 
correlated with the density of molecular hydrogen. This is 
the reason f or including the factor / c ,h 2 m Eq. [16] On the 
other hand, Idover fc Clarkl l|201lf l questioned a causal re- 
lationship between the star formation rate and the molec- 
ular hydrogen fraction. They argue that the observed cor- 
relation results form the necessity of effective shielding of 
star-forming regions from the interstellar radiation field. But 
this is in essence the effect that KMT09 describe with their 
model. For this reason, we also include the molecular hy- 
drogen fraction as a coefficient in the expression for the star 
formation rate. 

KMT09 implicitly account for the turbulent energy 
by assuming that molecular clouds are virialized. In addi- 
tion, the molecular cloud mass is determined by setting the 



Toomre stability parameter equal to unity. In Sect. [2721 we 
determine the mean cold-gas density p c , P a from an effective 
pressure balance, and we introduce a characteristic scale l c 
that is given by the thermal Jeans mass for this density. 
Since turbulence in the cold phase is generally supersonic, 
the local density of the gas greatly fluctuates. Therefore, 
we consider a statistical ensemble of overdense structures 
on the length scale l c . For convenience, we call these struc- 
tures clumps. The greater the overdensity relative to p c . P a, 
the smaller the critical density for gravitational collapse. 
For a given statistical distribution of density fluctuations, 
which we assume to be log-normal, the dimensionless star 
formation rate SFR Ci ff then can be calculated as proposed 
by PN11. 

3.1 Star formation efficiency 

By assuming a one-dimensional root mean square (rms) tur- 
bulent velocity dispersion a c , the virial p arameter of a clump 
at the mean density p c , pa is given by iBertoldi fc McKed 
( 1992) 
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7rGp c , P a V ' 



tc.ff 
^c.dyn 



(18) 



where t Cj dyn = 'c/(v / 3°"c) is the dynamical time scale. 
Since the statistical ensemble of clumps has to encompass 
the whole cold phase, not only regions of weak turbulent 
support, the rms velocity dispersion is given by the turbu- 
lent energy scaled down to the length scale l c . Hence, by 
substituting the scaling law (Eq. I12[) . the virial parameter 
can be expressed as 

10 e t (lg) 



«vir = 



TTGpcpa I 



2(1-17) 



/2<J 



PN11 argue that the overdensity in compressed shock 
layers is proportional to the square of the Mach number, 
M 2 (see Sect. 12.21 for the definition of M c )- By applying 
the Jeans criterion for the gravitational collapse of a com- 
pressed region within the cold gas phase, it follows that the 
critical overdensity ratio x CI a = Pc,crit/pc, P a is proportional 

tO Q v irA^c : 



Xorit = 0.0371a vir A4c 
0.0742 



(20) 



G7(7 - l)Pc, P a«c Z^-^V 



The constant of proportionality in the above equation 
is fixed by the definition of the Bonnor-Ebert radius (see 
PN11). For MHD turbulence, PN11 show that the x crit dif- 
fers by a factor /3 that specifies the ratio of thermal to mag- 
netic pressures. 

Since u c — const., the variation of the critical density 
Pcrit = 2; C ritPc, P a is solely determined by the second factor 
in Eq. I|20[> . For the two limiting cases of Kolmogorov and 
Burgers scaling, we obtain 

' e 2 z -2/3 r 4/3 if ,7 = 1/3, 

Pcrit oc <; 



if j) = 1/2. 



If the warm phase dominates (Z c <C I), then r\ w 1/3 be- 
cause, averaged over a region of size I, turbulence is mostly 
subsonic. From the scaling law e t oc l 2 ^ , it follows that 



Pcrit OC I 



-2/3 



The scaling behaviour of the critical density 
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follows from the steeper decrease of self-gravity with the 
clump size relative to the lower turbulent energy on smaller 
length scales. The assumption of Kolmogorov scaling is not 
at odds with supersonic turbulence within the clumps, be- 
cause supersonic scaling applies to length scales / < l c only. 
However, the assumption of a uniform velocity dispersion 
among both phases might break down for large A4 C - On 
the other hand, if the cold gas phase fills most of the vol- 
ume V — I 3 , r\ assumes a value greater than 1/3, depend- 
ing on M c - In the limit of high turbulent Mach numbers, 
Burgers scaling (et oc I) implies that p cr it becomes nearly 
scale-invariant □ In this case, however, a potential problem is 
that the turbulent pressure within the overdense cores (i. e., 
on the length scale of the shock-compressed layer, which 
is small compared to l c ) can exceed the thermal pressure. 
Consequently, the model overestimates the the mass that 
can form star in the limit of strongly supersonic clumps of 
size l c ~ I. Applying the model in numerical simulations, it 
has to be ensured that this case is sufficiently rare. 

As in KM05, the mass fraction per free fall time that is 
converted into stars is derived from the formula 



SFRc.ff = 



xp(x) dx, 



(21) 



where p(x) is the probability density function (pdf) of 
the mass density, and x = p c .ioc//5c, P a is the ratio of the local 
and mean densities in the cold phase. 

For isothermal gas, the probability density function 
(pdf) of the gas density is approximately log-normal (e.g.. 



(par) 01 the gas density is approximately 
IKritsuk et alJl2007l : IFederrath et alJl200a ) : 



p(x) dx - 



(2trt 2 ) 1/2 



exp 



(\n(x) - (ln(x))) 2 



2a 2 



dx, (22) 



where a 2 = — (ln(a;)} is the standard deviation of loga- 
rithmic overdensity. Log-normal fits to the density pdfs from 
the numerical simulations suggest the following empirical re- 
lation between a and the sonic Mach number: 



ln(l + b 2 M 2 c ) 



(23) 

As shown by IFederrath et all (|2010r i. the parameter b 
depends on the mixture of solenoidal and compressive forc- 
ing modes, which is specified by the weighing parameter £ 
of the Helmholtz decomposition of the forcing modes: 



, 1 2 
6 =3 + 3 



(i-C) 2 

1-2C + 3< 2 



(24) 



For solenoidal (divergence- free) forcing, f = 1. On the 
other hand, £ = for compressive (rotation-free) forcing. 
Substituting the log- normal pdf (|2"2"1) into Eq. (J2U, the di- 
mensionless star formation rate is given by 



SFR C 



1 + icrf 

2 2 



2 In (x c 



2 3/2 



(25) 



Numerical simulations of self-grav i tating turbulence 



(e.g. lKlessenll200ll: IFederrath et al.ll2008l : ICho fc Kimll201ll : 
IKritsuk et alj 12011 ) show changes of the high-density tail 
of the pdf, which affect SFR Cj ff .They find a power-law tail, 



2 In this case, the coefficient following from the assumption of 
spherical clumps would not be appropriate, but the scaling re- 
mains unaffected. 



wh ich is associated with self - gravit ating cores. Simulations 
by iBallesteros-Paredes et all (|201ll ) suggest that in a star 
forming cloud the pdf only develops a powerlaw tail at high 
densities over periods of > lOMyr, while the contribution 
of self-gravitating cores to the pdf is negligible in the ear- 
lier phase and, thus, the shape is close to log-normal. Since 
the model of PN11 for SFR Cj h is conceptually based on 
the turbulence-dominated phase, it is consistent to assume 
a log-normal pdf. An advanced formulation of the model 
might also account for the later power-law phase, but this 
would also require substantial modifications in the ansatz 
for SFR c> fj. We do not consider this in the present work. 

Furthermore, we assume a distribution of clump over- 
densities that is determined by the global rms turbulent en- 
ergy to estimate the fraction of collapsing gas in our PN11- 
like calculation of the star formation efficiency. This amounts 
to a separation of the density and velocity fluctuations. 
Strictly, the fraction of cold gas that can collapse should be 
calculated from the distributions of both the density and the 
turbu lent velocity fluctuations. As iHennebelle fc Chabrierl 
(2008) have already pointed out, however, this is far from 
trivial, and we do not attempt to solve this problem here. 



3.2 Molecular hydrogen fraction 

The formation of HVmolecules as well as their radiative de- 
struction are mostly heating processes. Because both rates 
are enhanced with density, overdense regions in a clump of 
cool but not molecular gas may be dispersed by this heating 
effect, before they possibly collapse gravitationally. So know- 
ing the fraction of molecular dominated gas in a clump, as a 
tracer for the fraction of gas that is not affected by effective 
radiation induced heating, is essential to correctly estimate 
the star formation rate. The fraction of molecular dominated 
gas in a cold clump is strongly dependent on shielding radi- 
ation, which dissociates HVmolecules easily, from its inner 
parts. Here we use a Stomgren-like approach similar to that 
iMcKee fc Krumholj l|2010l ) used. In low metallicity environ- 
ments this approach may lead to too high molecular fraction 
estimates, as reaction rates are too slo w to establish dissoci- 
ation equilibrium on short time scales (|Krumholz fc Gnedinl 
2011). But for our purpose this is fair enough, as we do 
not intend to track a whole chemical network of several 
sp ecies. Moreove r a sim ple chemical network model, like that 
of iGnedin et alj (|2009T h may have weaknesses, as it particu- 
larly in the case of large H2-fractions, which is of particular 
interest when lookin g at star formation , over estimates fur- 
ther H2-production (|Milosavlievic et al.ll20lih . Apart from 
that, this approach is not designed to compute the total 
fraction of molecular gas but the fraction that is molecular 
dominated, as we totally neglect molecular hydrogen in radi- 
ation dominated areas. Nevertheless we compare the results 
of this appoach to obervations in Sect. 16.41 
Assuming spherical clouds with diameter l c , one needs to cal- 
culate the radius Z c ,h 2 > a ^ which the production rate i?H 2 ,prod 
of H2 becomes greater than its destruction rate i?H 2 ,diss- 

The molecular fraction of cold gas then can be expressed 
as the ratio of the molecular volume in a clump oc ^c,h 2 an( ^ 
the total volume of the clump oc (l c /2) 3 : 



/c,H 2 = 



2/c 



(26) 
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where l c ,n 2 meets the condition 

-Z^H 2 ,prod 



^H 2 ,diss(rfc,H 2 



= 1, 



(27) 



where d c ,n 2 is the effective shielding layer thickness d, 
at a position inside the clump, where equation (|27[) is true. 
The iZVproduction and -destruction rates depending on d, 
assuming extinction of dissociating radiation of the outer 
regions of the cold clump is not sufficient, are given by (see 
iGnedin et aHl2009h 



i?H 2 ,prod — P°-P a (J p X ( — — rH 2 ,p,s + -^nj 2 ,p,g 

and 

-RH 2 ,diss(d) = IuSdustSll 2 — — ru 2 ,d 

m H 



(28) 



(29) 



respectively, w here C p = e CT is the clumping factor 
l|Gnedin et alj|2009l ) with a as defined in equation (|23[1 , Zq 
the solar metal fraction, 7"h 2 , p , s the H2-formation rate on 
dust surfaces, ?"H 2 ,p,g the Information rate in the gaseous 
phase, »"H 2 ,d the radiative dissociation rate, I v the intensity 
of the hom ogeneou s isotr opic dissociating radiation field rel- 
ative to the lDrainel (|l978l ')-fielcl. Sdust an d 5*h, are the shield- 
ing fa ctor s due to dust and H2 itsel f fsee lGlover fc Mac Low! 
(2007i) or lDraine fc Bertoldil (|l996l )): 



(30) 



Sdust = BXp ( — <7dust-^-Pc,pad 



1 - ^H 2 
>5H 2 — 7^— To + 



^H 2 



(i + x) 2 (i + x y. 



■ exp 



(-os a (l + x)*) (31) 



where 



X = /c,H 2 ,0Pc,pad/(mHK) 

with k = 5 • 10 14 cm~ 2 and 
/c,h 2 ,o = max (/c,H; 



(32) 



, -RH 2 ,prod/-RH 2 ,diss(rf = 0) J (33) 

( fc, Ho ,min « 10" 5 is the minimum molecular fraction 
in radiation dominated regions of the cold phase). 
In the centre of the spherical clump of diameter Z c the 
shielding layer has the same thickness for all directions, i.e. 
d = Zc/2. So if 

-Rh 2 , prod /Rn 2 , diss (d = Z c /2) 1 (34) 

holds, there is no molecular core in the clump, and thus 
/c,h 2 = 0. Otherwise there is one, which then is assumed 
to effectively block all dissociating radiation, trying to pass 
it. For a given position of scaled distance A = 21 /l c from 
the centre outside the molecular core (l c ,n 2 < I) the scaled 
effective absorption layer thickness 5 — 2d/l c is given by the 
mean of the absorption layer thicknesses over all sky O, but 
the solid angle of the molecular core 5 

S{\) = f S'(Q,)g(Q)dQ / [ g(Q)dQ , (35) 

Jo\S(X) I JO\S(X) 

weighted by the fraction of transmitted radiation, which 
is approximated by g = e~ s . The number of photons, that 
can possibly reach that position, is due to the cores shadow 
reduced by 

h. 

4-k 



(36) 



At the edge of the core A c ,h 2 = 2Z Ci h 2 /Zc half the sky 
is obscured (i.e. /„, shadow (A c ,h 2 ) = Iu/2). After integrat- 
ing/substituting out all angular dependencies we have 



•5max / rSn 



WA c ,h 3 ) = J s 5'g{8')A8' I g(S')dS' , (37) 



with Sxxxin = 1 — A c ,h 2 , 5 max = Jl - A J M and 

' 1 + A C 2 . H , - <5 



g(6') = 4itS' 2 e~ 6 1 



c,H 2 
'2 \ 2 \ ~ 2 



2A C 



(38) 



If the equations (|28[) to (|33l) are substituted into ()2T|) 
and using /^shadow (A c ,h 2 ) instead of J„, one obtains a tran- 
scendent equation for z = (x + 1) 5 : 

C.Ct^to + ^A-^- 1 ' (39) 



where 
C = 

D 



Pc^ f zr u 

lu, shadow rn 2 ,d m H V Z © n 2,P,t 



E : 



/h 2 ,o% ' 

2ftmy 
^cPc,pa/H 2 ,0 



(40) 



Eqn. (|39[) has a single solution for every given C, but 
only solutions in the range of z G [1 . . . jz max [ are allowed, as 
Z c ,h 2 would be greater than Z c /2 if z < 1 anc0 Z c ,h 2 ^ if 
z ^ max = (l + E- r )i. 

As 5 c ,h 2 (Ac,h 2 ) is bijective for A c ,h 2 G [0. . . 1], we can use 
its inverse A c ,h 2 (<5c,h 2 ) to compute the molecular fraction 



{0 if eq. flaUl true, 

(A C>H2 (1 - [z 2 - l]E)f if 1< z, 
1 else. 



(41) 



4 EVOLUTIONARY EQUATIONS 

4.1 Exchange of mass between the phases 

The effective growth rate of the stellar mass density is given 
by 



Ps.eff = pa — Ps,fb, 



(42) 



S(A) 



where the star formation rate p s is defined in Sect. |3j 
and p s ,fb is the rate at which gas is returned to the warm 
phase via core collapse supernovae (SNe II). 

In our model, p Sf fb is determined by a convolution of 
the past star formation rate p s (t — t') and the stellar initial 
mass function (IMF) d-/V*/dm* times the initial stellar mass 
ra*: 

f tc . , /x 1 dA» dm, , , . . 

7 p s (t — t) ~Wl dmT ~~d~F~ ' ^ 

where m* = m*(t', Z) is the initial mass of a star that 
explodes as a supernova after a lifetime t'. The integration 

3 Note, that the following case is already covered by an even more 
restrictive condition given in equation 11341 . 
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boundaries % to t e correspond to the lifetimes of 4OM0- 
and 8M0-stars, respectively. The IMF is normalised by the 
mean initial mass per solar mass 



i r° AN, ^ 

M* = -r—— I m*- dm*. 

M Q J dm* 



(44) 



For the functio n m*(t',Z), we use a parametrization 
(IRaiteri et al. 1996) of th e results comp uted by the Padova 
group ( Alongi et alJI 19931 ; iBressan et al1ll993h iBertelli et al 
119941 ). Furthermore, we assume the IMF bv lChabrierl (|200 



dA* 
dm* 



oc < 







for 0.1M Q < m, < 1M Q 
for 1M < m* < 125M Q 
otherwise, 



(45) 

where a,, c = 0.69 and m,, = 0.08 Mq. 
The fraction of heavy elements in the gas increases due to 
SN feedback. Assuming that the metal species in the ejecta 
have solar relative abundances, and, that the mass fraction 
of newly build up metals in the ejecta of a stellar population 
is independent of its initial metallicity at £ m as 0.1, we write 



d(Z P ) 
dt 



. . /~n 1 dA* dm, , , 



(46) 



with i = t — t ■ 

The rate of change of the fractional density of the cold phase, 
p c , is determined by the processes that are described in 
the following. The fir s t thre e processes are modelled as in 
ISpringel fc Hernquistl (|2003l ). For a schematic overview, see 
Figure [T] 



(i) Star formation reduces the mass of cold gas: 



dp c 
dt 



(47) 



(ii) Hot SN bubbles can evaporate cold clumps. Effec- 
tively, the energy that is injected by blast waves into the 
interstellar gas is instantaneously dissipated into heat on 
length scales that are much smaller than I. Since we cannot 
resolve the mixing processes, the dissipative heating and the 
heat conduction on these scales, we account for these pro- 
cesses by an evaporation rate of the cold gas, 



dp c 
dt 



Ap s ,tb, 



(48) 



where p B ,fb is defined by Eq . ()43l) . Following the analytical 
model of iMcKee fc Ostrikerl (|l977f ) for SN blast waves, the 
evaporation efficiency parameter A is given by 



P_ 

Pw,0 



lc 
lc,0 



>= >"(^-l l~) (^} ■ N!.l 



where p v = p — p c , and the length scale l c and the volume 
V c of the cold clumps are defined by Eqs. and ([9]). To ex- 
press variables in dimension-free form, we use the following 



scales: 

T = Tti = 1.5 x 10 4 K, 
uo = fcsTo/ '//mfl(7 - 1), 
Ao = usn/2mo, 
po = (itoh x 10.0 cm~ 3 , 

Pc,o = po, 
Pw,o = po — p c , o, 

lc,0 = Aj jC (lt c , p c ,pa,o)- 

(iii) The cold phase gains mass from the warm phase via 
radiative cooling if the gas is thermally unstable: 



dp c 
dt 



(1 — ,/th)/9wA e ff 



(50) 



where the effective cooling rate A c ff is defined below in 
Sect. 14.21 and / t h is the thermal stability indicator that 
switches on/off terms in the model equations that are re- 
lated to the thermal instability: 



/th := 



if conditions 
else. 



a)-(c) are satisfied, 



(51) 



The warm neutral gas is treated to be thermally unstable, 
in the meaning of currently separating into two pases due to 
cooling, if following conditions are met: 

(a) The net effect of radiative cooling, Lyman contin- 
uum radiation, UV background radiation, and turbulent 
dissipative heating must decrease of the thermal energy, 
i.e. A cff > 0. 

(b) The warm gas density p w , P a must exceed 
O.l^tmHcm -3 , since this is roughly the minimum density 
for thermal instability according to the equilibrium cool- 
ing curve. 

(c) Furthermore the gas must be largely neutral to be 
thermally unstable, thus we assume an upper temperature 
threshold Tti = 1-5 • 10 4 K for the cooling instability of 
the warm gas. 

(iv) Massive stars of spectral class O and B are strong 
emitters of radiation in the far ultraviolet. In particular the 
photons in range of the Lyman continuum (Lyc) deposit a 
significant amount of energy iLyc per photon (See Sect. 16.3) 
in the gas, as they are absorbed and then reemitted as Ly- 
man q (Lya?) photons. The number of Lyc photons emitted 
by young, massive stars per unit volume and per unit time, 
ALyc,ioc, is computed from a convolution of the past star for- 
mation rate p s and the specific emission rate of Lyc-photons 



iVLy C ,ioc(t) = / Ps (i')™Lyc (*-t'i Z)dt' , 
Jo 



(52) 



For nLyc(£ — t',Z), we use an analytic fit to data from 
evolutionary synthe sis models of a simple star population 
l|Kotulla et al.ll2009l) . Some fraction /i oa k of these photons 
may leak into the environment of the reference volume, while 
ALyc,ext photons from external sources may get in. The effec- 
tive number of Lyc photons per unit time and unit volume 
that actually ionize hydrogen is then given by 

A L yc = (A L yc,loc + ALyc,OXt)(l ~ /leak), (53) 

where /i ea k — exp(— anXpl/mYi) with the ionization 
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crossection <th — 6.3 • 10 cm . In case of average hydro- 
gen number densities pX/mu ^ 1 cm -3 and length scales I 
of at least a few parsec, /i ea k is negligible small. It is likely 
that the reference volume is surrounded by an environment 
of comparable density and size, if sources of Lyc radiation 
are located there, thus iVLyc.ext is negligible, too. Hence, we 
set A^Lyc = A^Lyc.ioc and assume that every hydrogen atom 
has the same chance to absorb a photon. The radiative evap- 
oration rate is given by: 



dp c 
dt 



(u„ — u c )p ' 



(54) 



(v) Turbulent energy is dissipated into thermal energy at 
a rate that is given by 



e = C* E 



3/2 
St 

I 1 



(55) 



where C c is about unity (see ISchmidt et al ] l2006t SF11). 
In our two-phase model, turbulent dissipation heats the gas 
in the warm phase, but we assume the temperature of the 
cold gas to be constant. To compensate turbulent dissipation 
in the cold phase, we simply transfer an equivalent amount 
of mass from the cold to the warm phase. Since turbulence is 
assumed to be homogeneous on length scales smaller than I, 
the amount of energy dissipated in the cold phase is m c edt 
over an infinitesimal time interval dt. Setting this equal to 
the increase of thermal energy (w w — u c )dm c if the mass dm c 
is transferred to the warm phase and substituting m c = p c V, 
we obtain 



dpc 
dt 



n 3/2 

C t e t > p c 

l(tlw — U c 



(56) 



(vi) If the cold gas forms small compact clumps embedded 
in the warm phase, i. e., l c <C I, we can model interactions 
between the clumps as collisions. Since collisions cause a 
certain mass loss of the cold phase by turbulent mixing and 
heating, we write 



dp c 
dt 



i3 

~£ccPc,pa^c, collie 



(57) 



The effect of clump collisions on the cold gas density is 
modelled by the efficiency parameter e cc and the the collision 
rate 

r V C ,C0ll 



^"c.coll 



n c V 



^c,fre 



(58) 



Setting the average volume of a cold clump equal to irl^ /6, 
the number density of the clumps is n c ~ {6V c /nl 3 c )/V and 



the mean free path Z Ci f roc = (tyI 



The rms velocity of 



the clump motion in the surrounding warm medium can be 
estimated from the square root of the turbulent energy et, 
corrected by the internal velocity dispersion a\ of the clumps 
(see Eq.[l2j): 



(2et - 3a, 



2\l/2 



2e t 1 



2,-1 



1/2 



With the above definitions, it follows that 



^"c.coll 



36K 2 

TVIW 



2e t 1 



2>! 



1/2 



(59) 



(60) 



where V c and l c are given by Eqs. © and (0 , respectively. 




Figure 1. Scheme of the exchange of mass. The mass budgets 
are depicted in yellow (p s ), red (p w ) and blue (p c ), where the 
molecular mass (/ c ,h 2 Pc) in darker blue resides inside the cold 
gas mass. Arrows illustrate processes transferring mass from one 
to another budget. 



Collecting all six contributions, the evolutionary equa- 
tion of the cold phase density can be written as 



p c = — p s — ^4ps,fb — 



-(1 — JthJPwAcff H 1 ; 



36e cc p CiPa V e 2 
irl c 



2e t 1 - KS 



2 n 



1/2 



(61) 



and the change of the gas density in the warm phase 
follows from mass conservation (p w + p c + Ps.eff = 0): 



pw 



-pc — Ps.cft'- 



(62) 



4.2 Exchange of energy between the phases 

In nu merical simulations of the r mally bistable turb ulence 
(e. g., lAudit fc HennebellellioTol ; ISeifried et alj|201ll ). most 
of the gas in the cold phase is close to the isothermal equilib- 
rium branch of the cooling curve. This is mainly caused by 
the higher opacity of the dense, cold gas that lowers the ef- 
ficiency of radiative cooling. Other processes that affect the 
gas temperature in the cold phase, such as the gravitational 
collapse of dense regions or chemical reactions, are not ex- 
plicitly treated in our model. Consequently, the specific ther- 
mal energy of the cold phase u c is assumed to be constant. 
We set the temperature to a fiducial value T c = 50 K, cor- 
responding to the lower cutoff of the cooling curve in galaxy 
simulations. 

The specific thermal energy of the warm phase, on the 
other hand, is changed by the processes that are discussed in 
Sect. 14. ll The effects of these processes on w w are as follows 
(see also Figure [2}. 
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(i) SNe heat the warm gas and transfer gas from the cold 
to the warm phase via evaporation. The energy release per 
unit mass is usn ~ 6 ■ 10 49 erg/A/e (see Sect. I|6,3j l), As- 
suming that a certain fraction esNWSN of the feedback is 
non-thermal (see Sect. I4.3|) . we have 



d(p„ 



dt 



[(1 — £Sn)msn + Au c ]/9 s ,fb. 



(63) 



SN 



The efficiency parameter A is defined by Eq. (j49J). The 
rate of change of the warm gas density due to SN feedback 
is given by the mass ejection from SNe and the evaporation 
of cold gas: 



dp w 
dt 



= (f + A)p B , {b , 



Combining the above equations, it follows that 



dii w 



[(1 - ESNWSN + Au c - (f + A)u w )] 



Ps.fb 



(64) 



(65) 



(ii) If the warm phase is thermally stable, the warm gas 
cools (or heats) at a rate given by its effective cooling func- 
tion A c ff. Once the thermal instability sets in, gas in the 
warm phase is converted into cold gas without changing 
the temperature of the remaining warm gas. We also as- 
sume that the cooling instability produces turbulent energy 
with an efficiency e t t relative to the cooling function. Con- 
sequently, the total change of the internal energy density of 
the warm phase can be written as 



d(pw^w) 
dt 



a i d ^ w 
= — p w A off + u c — — 

dt 



(66) 



-(1 



/th)ettPwA e ff. 



Since the rate of change of the warm gas density due to 
cooling is given by Eq. (|50[1 multiplied by minus one, we 
obtain 



du„ 
~dT 



= — /th-A-efl — (1 — /thJCttAeff, 



The effective cooling rate A e g is defined by 



A c ff = A fl 



Lyc 



(67) 



(68) 



where A ra d is the specific radiative cooling rate. In this 
model, we use a tabled atomic cooling function, computed 
using the photo-i onisation package C loudy (version 08.00), 
last described bv lFerland et al.l (|f 998l l. without considering 
molecules or dust. Tpah is the photo-electric heating rate 
l|Wolfire et al.lll995l ) due to the external radiation field I v 
modified by a factor of Z/Zq, and e is the turbulent dissi- 
pation rate per unit mass ()55[) . The volume rate of heating 
by Lyc photons is given by 



d(p„ 



dt 



AT I d Pw 

JVLyc^Lyc + «c "~ ~jT~" 



Hence, the specific heating rate is 



r Lyc — 



du w 
~dT 



(69) 



(70) 



(iii) Since cold gas is transferred to the warm phase by 
clump collisions, we have 



d(p„ 



df 



dpc 
dt 



(71) 



The corresponding rate of change of the specific energy is 
given by (see Eq. 



du w 
~~dT 



^CC (^V, 



where e cc is the efficiency parameter of the collisions, and 
the collision rate r CiCO n is defined by Eq. (|60p . 

Adding up the contributions (i) to (v), the dynamical 
equation for the thermal energy of the warm phase becomes 

1 Ps.fb 



= [(1 — esn)«sn + Au c — (1 + A)w w )l — — 

Pw 



[/th + (1 - /th)ett]A e ff 

/ \ Pc,pa^*c,coll^c 



(73) 



4.3 Turbulent energy production and dissipation 

To formulate an equation for the turbulent energy, we as- 
sume that energy is injected at constant rate E by large-scale 
forcing. The rate of energy injection determines the flux of 
kinetic energy that is transported through the turbulent cas- 
cade from larger to smaller scales. For purely hydrodynamic 
isotropic turbulence, the energy flux is independent of the 
length scale and equal to the dissipation rate e in statisti- 
cal equilibrium. Applying the method of (adaptively refined) 
large eddy simulations, E can be computed from the SFf I 
closure for the compressible turbulent energy cascade. For 
the one-zone formulation of our model, we simply express 
E in terms of the typical magnitude of the turbulent veloc- 
ity fluctuations V induced by the turbulent cascade on the 
length scale I: 



E = C e p 



V 3/2 



(74) 



For pure hydrodynamical turbulence, E = e and e t = 
0.5V 2 in equilibrium. 

Neglecting turbulent diffusion and collecting the terms 
that exchange energy between the gas phases and turbu- 
lence, the following rate equation for the turbulent energy 
results: 



e t = (esNWSN ~ e *) 1" I 1 _ Ah)£ttA c ff — 

P P 



+ C E -E- 

p I 



3/2 



(75) 



The three sources of turbulent energy production are 
SN feedback on length scales comparable to the cooling 
instability and the turbulent energy cascade. The two ef- 
ficiency parameters esn and ett determine the coupling of 
the unresolved processes to the turbulent energy. In numer- 
ical simulations, in which / corresponds to the grid scale, 
these parameters have to be chosen appropriately. The cru- 
cial problem is that, in contrast to the cascade of turbulent 
eddies in the inertial sub-range of isotropic turbulence, SN 
feedback and the cooling instability are not self-similar. We 
can only assume that particular efficiency parameters apply 
to certain ranges of scales. One option is to use small-scale 
simulations of the interaction of SNe blast waves with the 
interstellar medium and thermally bistable flows in periodic 
boxes to estimate these parameters. On the other hand, the 
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Figure 2. Scheme of the exchange of energy. The ellipse depicts 
the energy content of the gas, which is separated into three bud- 
gets: the thermal energy of the warm gas p w ^ w (red) and the 
cold gas p c u c (blue) and the non-thermal, turbulent energy pet 
(green). Cold phase thermal energy can only change by loss or 
gain of mass, as u c = const., this is shown by dashed arrows, while 
processes not intermixing cold gas with any phase have solid ar- 
rows, star formation removes mass from the cold gas, along with 
its energies, which is shown in dot-dashed arrows. 



model can be calibrated a posteriori in large-scale simula- 
tions such that observational constraints are met. We also 
show in Sect. [6l that one-zone calculations can be utilised 
to find reasonable choices of the efficiency parameters. 

Combining Eqs. (|61I62I73I75[) and u c = const., the fol- 
lowing law of energy conservation equation is obtained: 

d(p(u + e t )) 



dt 



= - p s (u c + e t ) + p s ,fbWSN 

— pwA rad + p w r P AH ( 76 ) 

+ A^Lyc^Lyc + S. 



5 LIMITING CASES 
5.1 One-phase medium 

If most of the gas cools down to temperatures close to T c , 
the separation between a dynamic warm phase and a star- 
forming cold phase breaks down. Assuming that the warm 
gas can also form stars if « w u w ,mm, where it w .mm u c , 
the star formation law (11511 becomes 



(77) 



where p c ,pa = P and p — — p s ,cff- In this article, our 
focus is on a statistical description. By choosing a length 
scale I 3> l c — ^j{p), the star formation time scale t s can 
be calculated as in Sect.0 When applying this model as an 
SGS model in numerical simulations, however, the Truelove 
criterion requires I ^ Zj(p)/4. Consequently, an alternative 
parametrization of the star formation time scale has to be 
applied in the one-phase limit. This is left for future work. 
Practically, this case will occur in simulations of individual 
galaxies with high resolution, but usually not in cosmological 
simulations, where cold-gas clumps are sufficiently below the 
resolution limit. 



Since the cooling instability vanishes, and the exchange 
of energy between the phases as well as the collisions of cold 
clumps become meaningless, the equation for the specific 
thermal energy becomes 



ii„ ~ [(1 — e SN )MSN 



Uc)] — 



(78) 



The subtraction of u c in the factor that is multiplied 
with the feedback rate results from the removal of the energy 
in the gas that forms stars from the gas energy budget (see 
Eq. I76|) . Practically, we can neglect the difference because 
usn S> u c . The thermal energy equation is complemented 
by the simplified turbulent energy equation: 

3/2 

(79) 



1 e 

et = - (eSNfiSNps.fb + E) — C e — — 

p I 



It is instructive to consider the asymptotic limit u w ~ 
u c = const. Then we can set u w ~ 0. For net heating (A e g < 
0), this equation cannot be fulfilled because both terms are 
positive. If A e fl > 0, on the other hand, the feedback rate is 
approximately given by a balance between thermal heating 
by SNe and turbulence production by cooling: 

Ps.fb ==: -n : • (80) 

(1 — esn)msn 

Since the p Sj fb is related to the star formation rate (|77|l 
via Eq. (|43[) . the above equation imposes a condition on the 
effective cooling rate so that u w — u c = const. 



5.2 Equilibrium solutions 

Of particular interest is the case of self-regulation, for which 
the star formation rate is low and nearly constant: p s ~ 
p s ,eq = const. A low star formation rate means that changes 
in the gas density are negligible in first-order approximation. 
In addition, we assume that the temperature of the warm 
phase and the specific turbulent energy are approximately 
constant in the self-regulated regime and that the cooling 
instability is active (/ = 0). For simplicity, we neglect clump 
collisions. Setting w w ~ and e t ~ in Eqs. ([73} and (|75[) , 
respectively, it follows that 

[(1 — £Sn)usn + Au c — (1 + A)Uw,cq)]^^ — £ttA cff ~ 0, 

Pw 



(81) 



3/2 



, \ Ps.fb . Pw £ „ e t,cq 
(eSN«SN — e+eq) V CttA c ff 1 C e — ; — ~ 0, 

p p p I 

(82) 

where u w ,cq and et.eq are the equilibrium values. 

Equation (]82fl imposes a condition on the feedback rate. 
By substituting the effective cooling rate (|68[) for A e fi, we 
obtain 



(esNttSN — et.cq)^^ — C E (p + ettPw)-^ 2 - _ s 
P I 

+ £ttPw(rpAH + TLyc — A ra d). 

(83) 

For any reasonable choice of parameters, et, eq <C 
esnWsn- Consequently, a solution exits only if 

3/2 

g 

C e (p + ettPw)— j— ^ S + ettPw(A rad — Tpah — TLyc). 



3/2 
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The right-hand side is always positive if the cooling insta- 
bility is active, because A ra( j — Tpah — IYyc ^ Aes > 0. For 
this reason, the above constraint implies that there must be 
a minimal turbulent energy for self-regulation: 



mm et,. 



[E + ettPw(Arad — TpAH — T 



Lye, 



C e (p+ ettPw) 



2/3 



This an important implication of the multi-phase model. 

The stellar feedback rate is also a constant in the self- 
regulated regime. This follows immediately from Eq. (|43[) 
for a constant star formation rate. Then 



where 



Ps.fb — Pps 



1 AN, dm, 
M»dm„ At' 



At'. 



(84) 



(85) 



In principle, one could invert the equations for a given 
star formation rate, p s , cq , and the effective pressure equi- 
librium between the phases, to obtain the cold and warm- 
gas densities. Due to the high non-linearity of these equa- 
tions, particularly the molecular hydrogen fraction, this is 
very difficult in practice. It is easier to search for equilib- 
rium solutions by computing the full set of rate equations 
and identifying solutions that are close to equilibrium values 
satisfying Eqs. (|81|l . (|82|) . and (|84|) in certain time intervals. 

By substituting Eq. (152)) into Eq. flST}, the following 
expression for the equilibrium energy of the warm phase 
can be obtained: 



MSN A 

+ — TUc 



1 



1 + j4 e t , cq 



(l + A)Pp\ c 



3/2 
t.cq 



(86) 



For E = (no turbulence feeding by instabilities on 
length scales greater than I) and e t , eq = (turbulent energy 
is neglected), the SH03 equilibrium solution for u w , cq results, 
with the exception of the factor A /(I + A). Since A » 1, 
however, this factor is very close to unity and the result is 
practically the same. In our model, the cooling instability 
produces turbulent energy on top of the turbulent cascade. 
Thus, e t ,eq > (E;/C £ p) 2/3 (see Eq.|ST) and 



Mw, eq < MSH : = 



+ Mc 



(87) 



As a consequence, we expect that the temperature of 
warm gas close to equilibrium decreases in numerical simu- 
lations with a turbulence SGS model, because a fraction of 
the energy is in non-thermal form. 



6 MODELLING THE EVOLUTION OF A 
SINGLE ZONE 

The set of six coupled nonlinear differential equations (|42l 
EDE3E3E51HS1), as defined in Sect.g] describe how the gas 
in the reference volume evolves with time. By numerically 
integrating the model equations over closed boxes (i. e., sin- 
gle zones), we obtain statistical models for a wide range of 
initial conditions and parameters. These models also allow 
us to find equilibrium states, for which the star formation 



rate is small and nearly const, (see Sect. 15. 2[) . To charac- 
terise the star formation rate, we define a dimensionless star 
formation efficiency by 



£ff 



Pstff 

9 



This is the fraction of the total gas mass in the ref- 
erence volume that is turned into stars over a free fall 
timescale tg — ^/37r/32Gp. It is to be distinguished from 
SFR c .ff, which is sometimes also called star formation effi- 
ciency. However, SFR c ,ff specifies the fraction of cold molec- 
ular gas converted per free fall time (see Sect. [3}. 

If not stated otherwise, the following standard model 
parameters are used: 

V = 1/3, 
e tt = 0.025, 
e cc = 0.0, 

I — 15 pc, 

We comment on the choice of these units in Sect. 16.31 
Furthermore, metal enrichment is turned off by setting £ m = 
0. We specify the strength of the turbulent energy injection 
E by relating the velocity scale V (see Eq. I74[) to the speed 
of sound at the temperature Tti = 1.5- 10 4 K, corresponding 
to the maximum thermal energy uti of thermally unstable 
gas. This results in the Mach number 



b 


= 2/3, 




= 0.085, 


^Lyc 


= 0.1 eV, 


floss 


= 0.6. 



Me = 



V_ 

CTI 



7(7 — 1)mti 



1/2 



SI 
Op 



1/3 



(89) 



as a basic parameter for the external (large-scale) energy 
injection relative to the maximal turbulence production by 
the thermal instability. We use the standard value My, — 
0.2. 

In the following sections 16.11 and 16.21 sample evolutions 
of gas are discussed. The gas, as well as the formed stars, 
in these evolutions are confined to the reference volume, i.e. 
nothing but energy enters or leaves the volume. If the ref- 
erence volume was embedded in a more realistic, inhomoge- 
nous, dynamic environment, like in an isolated disk galaxy 
simulation, it would exchange mass with its environment, 
as gas pressure is subject to strong variation, depending on 
the enviroment this may lead to substantial in- or outflows, 
and stars may leave the area due to their drift. This is be- 
yond the scope of a single zone model, but to demonstrate 
the effects of basic model parameters it is still useful. For 
the equilibrium solutions discussed in Sect. 16.31 however the 
latter is no objection. 



6.1 Feedback sequence 



As described in Sects. 14.11 and 14.21 stellar feedback is de- 
termined by the evolution of the stellar populations within 
the reference volume (e.g., SN feedback starts when the first 
SNe II light up). As an example, we consider the evolution 
of gas with number density n = p/(pmu) = 75 cm -3 , solar 
metallicity, and no turbulent energy injection (i.e., E = 0). 
To single out the effects of the stellar feedback, we artificially 
suppress star formation after one Myr has passed. The re- 
sults are shown in Fig. [3] The evolution starts in an equilib- 
rium between heating and cooling without star formation. 
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Figure 3. Demonstration of the feedback sequence: Star for- 
mation is only active for t £ [0 ... 1] Myr (period between the 
dotted vertical lines) and set to zero at later time. The upper 
panel shows the mass fractions of the different phases as area fill- 
ings (stellar fraction p s /p, contribution of the young Lyc-emitting 
population p s ,Lyc/Pi warm fraction p w /p, neutral cold fraction 
p c (l — / c ,H 2 )/p> an d molecular fraction / c ,H 2 Pc/p). m the lower 
panel, the evolution of the star formation efficiency eg, the rms 
Mach number of motions in the cold phase M c , the size of cold 
clumps l c , the thermal energy of the warm phase u w , and the 
turbulent energy et are plotted. 

Then star formation is activated within the interval indi- 
cated by the vertical dotted lines. Lyc-heating by the newly 
formed stars quickly lowers the star formation efficiency eg 
defined by Eq. (|88|l . as the the cold phase gas fraction de- 
creases, Z c drops while / c ,h 2 grows, and the warm-gas pres- 
sure rises. But the specific energy tt w remains low because 
radiative cooling dominates in the warm phase. When Lyc- 
heating begins to decrease (vertical long-dashed line), more 
warm gas can cool down and / c ,h 2 drops due to the shift in 
the pressure equilibrium between the phases. After the start 
of SN feedback (short dot-dashed vertical line) u w and p w 
rise quickly, and the ensuing compression of the cold phase 
results in an increase of /c,h 2 - The kink in the evolution of 
l c and A4 C that can be seen shortly after the onset of feed- 
back is caused by the production of turbulent energy, which 
is slightly lagging behind. As SN feedback becomes weaker, 
u w , e t and / c ,h 2 are declining. Once all SN- progenitors are 
gone, SN feedback ceases (long dot dashed vertical line). 
Then the warm gas cools down to the initial temperature 
(reaching it at the short dashed line), and all quantities and 
gas fractions evolve back to the values for the equilibrium 
of heating and cooling without star formation. 

6.2 Dynamical evolution 

The main parameters of our model are the initial total gas 
density n, the metallicity Z and the forcing Mach number 
A^e, corresponding to the turbulent energy injection rate 
E. In this section, we describe the dependence of the phase 
evolution and star formation history on these parameters. 
In addition to p, Z, and A4e, several coefficients determine 
the relative contributions from unresolved processes. The 
influence of these coefficients is discussed in Sect. 16.31 



Three sample evolutions are plotted in Fig.|4]for the ini- 
tial number densities n = 35, 50 and 75 cm -3 , Z = Z Q , and 
AIe = 0.2. For the lowest density, no star formation occurs 
at all, except for a negligible fraction in the very beginning 
(the spurious SN feedback produced by these stars causes 
the kink in the cooling curve of u w at t ~ 4.2 Myr). The 
final /9c/pw is determined by the equilibrium between turbu- 
lent dissipation, photoelectric heating and radiative cooling. 
The final et is fixed by the equilibrium between the pro- 
duction of turbulent energy by the thermal instability and 
large-scale injection and turbulent dissipation. For higher 
values of n (middle and right panel of Fig. 3]), on the other 
hand, a markedly different evolution can be seen. Following 
an initial transient phase that ends after about 10 Myr, a 
stationary mode of star formation is entered, in which the 
star formation efficiency eg is a few per mil. During the 
transient phase, there are three more or less distinct max- 
ima of eg. This behaviour can be understood as follows. 
The initial rise of p w /p due to Lyc-heating by the first stars 
causes a compression of the cold gas. This results in an in- 
creasing molecular hydrogen fraction and, thus, an enhance- 
ment of eg. Depletion of the cold phase reverses this trend. 
As Lyc-heating fades out, the star formation efficiency in- 
creases again, resulting in a second, but weaker peak. Due 
to the cooling of the warm phase, which exerts pressure on 
the cold phase, p c , P a is declining and eg is lowered. The first 
SNe raise w w and the subsequent growth of turbulent energy 
causes the third maximum. 

By comparing Figs. [4] and [5] where the latter fig- 
ure shows plots for three one-zone models with an initial 
density n = 75cm -3 , but different metallicities Z/Zq — 
1.0, 0.8, 0.75, one can see that lowering Z has effects on the 
evolution similar to those of lowering n. This is simply a 
consequence of the dependence of the following processes on 
the density of metals, which is proportional to nZ: 

• A ra d is primarily determined by metal line cooling. 

• The major contribution to the Lb-production rate (see 
Eq. I28[l is the formation of molecules on dust grains, whose 
abundance is assumed to be proportional to Z. 

• The absorption of Lb-dissociating radiation outside of 
molecular cores is dominated by dust extinction. 

Consequently, / c ,h 2 and the relative mass fractions of phases 
are sensitive to Z, while quantities related to gravity are only 
indirectly affected. For this reason, l c , tg and SFR c ,ff remain 
almost unaffected when varying Z. 

Next, we consider the influence of turbulence driving. 
Generally, a higher production rate E increases the turbulent 
energy et and damps the peaks of star formation during the 
initial transient phase. Consequently, the stationary phase is 
entered earlier and more smoothly, as one can see in Fig. [6] 
Even in the absence of turbulent energy injection {Mt, = 0), 
the turbulence generated by the cooling instability and SNe 
plays an important role in limiting the star formation rate 
(see the left panel in Fig. [6} . This is caused by a decrease 
of p c , P a, as indicated by the growth of l c . Small A4e do 
not cause a major increase of e t , M c or l c (middle panel 
of Fig. [6]) compared to Ms ~ 0, although u w settles at 
a slightly higher level due to the increased energy input. 
However, if the production of turbulence is dominated by E, 
the turbulent contribution to the effective pressure becomes 
comparable to the thermal pressure in the warm phase (et — 
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Figure 4. One-zone models for different initial gas densities at solar metallicity an moderate external driving Ms = 0.2. See the legend 
in Fig. [3] for a definition of the plotted quantities. 
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Figure 6. One-zone models for different rates of turbulent energy injection at fixed initial density n = 50 cm 3 and solar metallicity. 
See the legend in Fig. [3] for a definition of the plotted quantities. 
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it w in the right panel of Fig. [6j. As a consequence, the cold- 
phase pressure increases relative to the warm-phase pressure 
(see Eq. [9j), and p c ,pa decreases. This results a lower star 
formation efficiency eg. Apart from this effect on eg for e t ~ 
u w , turbulence also affects the molecular fraction / c ,h 2 and 
the star formation efficiency in the molecular gas, SFR Ci g. 
These sub-dominant effects are discussed in Sect. 16.31 

6.3 Equilibrium star formation efficiency 

Although the transient phases discussed in Sect. 16.21 shed 
light on the complex interplay between the various physi- 
cal processes contributing to the multi-phase dynamics, the 
approximately statistically stationary regimes of star for- 
mation are most relevant for applications. To numerically 
determine equilibrium solutions, we assume a small star for- 
mation rate of 1 % as initial condition and integrate the rate 
equations until the relative temporal variance of p s over a 
feedback period (>40 Myr) becomes less than 10 -4 or the 
star formation rate approaches zero (in the latter case the 
solution is not considered to be an equilibrium solution). We 
also check if the equilibrium conditions (|81[) and (|82[) are ful- 
filled with a relative accuracy better than 10 -4 . In contrast 
to the dynamical evolutions discussed before, we keep the 
total gas mass constant in the course of the integration, be- 
cause we want to obtain equilibrium solutions for fixed gas 
densities that do not depend on the gas consumption by 
star formation during transient phases prior to the statisti- 
cally stationary states. In numerical simulations, where the 
model describes a certain grid zone, one can think of the gas 
being replenished by neighbouring regions. Of course, such 
conditions will be met only to a certain degree and for a lim- 
ited period of time. Nevertheless, the equilibrium solutions 
are useful to understand the behaviour of the system, and 
these solutions can be utilised as approximations to the star 
formation under quasi-stationary conditions. 

Gas density and metallicity 

Many star formation recipes used in astrophysical simula- 
tions assume, inspired by observational Kennicutt-Schmidt- 
relations, that stars are formed with a fixed efficiency per 
free fall time if a certain density threshold is exceeded. Our 
model shows such a behaviour under the condition of low 
turbulence driving E, as the star formation efficiency es sat- 
urates quickly above a metallicity-dependent density thresh- 
old (see Fig. \7\ left panel). The dependence of es on n is 
governed by the molecular fraction / c ,h 2 , while SFR Ci ft re- 
mains approximately constant for all n 2> 1 cm -3 . A similar 
dependence on the metallicity can be seen in Fig. [7] right 
panel. This is a consequence of cooling, tb-production and 
extinction of Hb-dissociating radiation being mainly depen- 
dent on nZ. 

For the zone evolutions with constant total mass (gas 
and stars) shown in the left panels of Figs. [4] and the star 
formation efficiency vanishes, although n and Z are above 
the threshold values for star formation following from Fig. [7) 
This can be understood as a consequence of the different ini- 
tial conditions. Due to the lack of stellar feedback, the gas 
phases are never pushed into a star-forming regime in the 
former calculations. In a numerical simulation, the heating 



of gas by nearby stars could trigger star formation in so far 
inactive gas. Otherwise the gas will remain cold and inac- 
tive as long as it does not become dense enough to form 
molecular cores and to start star formation. 

External turbulence driving 

Low turbulence intensity maintained by external driving 
(Ms < 0.2) does not change the shape of eg(n) apprecia- 
bly. Fig. [8] (left panel) shows that moderate values of Ms 
increase the density threshold and lower the saturation level 
of eg slightly. As Ms approaches unity, however, star for- 
mation becomes more and more inhibited. Turbulence also 
smears out the density threshold. In Fig. [8] (right panel), four 
regimes can be identified, in which the effects of turbulence 
injection differ with increasing Ms- 

(i) For small Ms, the additional heating of the gas by 
turbulent dissipation partially counters the effects of turbu- 
lence suppressing star formation that become dominant for 
stronger Ms- 

(ii) For higher Ms, the production of e t is dominated 
by the turbulent cascade and e% roughly follows M\. As a 
consequence, turbulent pressure contributes significantly to 
the pressure balance between the phases, and p c , pa decreases 
(/9c, P a — > p in the limit of large et). The lowering of the cold- 
gas density results in the steep reduction of eg as Ms rises. 

(iii) For stronger turbulence intensity, p c , P a is low but the 
turbulent broadening of the density pdf becomes important 
(enhancement of HVproduction by the clumping factor C p 
and dependence of SFR c ,fj on the pdf; see Sect. I3.l| l. This 
effect can clearly be seen for initial densities lower than 
50 cm -3 , where a second maximum of es can be discerned. 
This is the regime, in which star formation critically depends 
on the properties of self-gravitating turbulence. 

(iv) If Ms increases further, the growth of of the mini- 
mum overdensity for star formation in the cold gaS, £crit OC 
e\ (see Eqns. and |2"5")) dominates, and es asymptotically 
falls off to zero. 

Because of the dependence of SFR Cj g and / c ,h 2 ° n the 
density pdf, the saturation level of es and the density thresh- 
old are significantly affected by the weight of compressive 
forcing modes relative to solenoidal modes. Assuming that 
the width of the density pdf is given by a « log (I + M c° 2 ), 
where M c is the the rms-Mach number of turbulent mo- 
tions in the cold gas and b varies between 1/3 for purely 
solenoidal forcing and I for purely compressive forcing, we 
obtain the equilibrium solutions plotted in Figs. [5] and QJj] 
(right panel). The nature of turbulence driving in the inter- 
stellar medium is still a matter of debate. Moreover, the 
mixture of solenoidal and compressive modes is likely to 
be scale-dependent. Here, we adopt the intermediate value 
b = 2/3 as default. 

As explained in Sect. 15.21 turbulence generally decreases 
the temperature of the warm gas in equilibrium (see Eg. 1861) . 
Fig- HU shows that this is roughly a 10 % effect for densities 
higher than the threshold for star formation. In a certain 
sense, this is the deviation from the SH03 equilibrium so- 
lution ush (see Eq. [87)) . However, the value of usn used by 
SH03 and the coefficient of the first term is different from 
our definition so that the actual difference is greater. As 
one can see in Fig. 111! there is a noticeable deviation even 
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Figure 7. Equilibrium star formation efficiency per free fall time, eg, versus number density n for different metallicities Z/Zq = 
1.0, 0.8, 0.6, 0.4, 0.2 from dark to light blue (left panel), and eg as function of Z/Z Q for different n = 200, 400, 75, 50, 25, 42.5 cm -3 
from dark to light blue (right panel). 
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Figure 8. Equilibrium star formation efficiency per free fall time, eg, versus number density n for different turbulent energy injection 
rates S (left panel) in terms of AIe = 0.0, 0.2, 0.4, 0.6, 0.8, 4.0, 4.2 from dark to light blue (see Eq. I89H . and eg as function of Ms f° r 
different number densities n = 27.25, 32.5, 35, 37.5, 50, 400 cm -3 from light to dark blue (right panel). 



0.008 
0.006 
$ 0.004 
0.002 
0.000 








50 100 
n [cm" 



150 



200 



0.008 
0.006 
'0.004 
0.002 
0.000 



0.4 0.5 0.6 0.7 
b 



0.9 



1.0 



Figure 9. Equilibrium star formation efficiency per free fall time, eg, versus number density n for different values of the compressive 
factor b in the cold phase (left panel). From light to dark blue, b = 4/3, 0.5, 2/3, 0.75, 4.0. eg over turbulence forcing parameter b (right 
panel) for the number densities n = 200, 400, 50, 37.5, 34.25, 25 cm -3 (from dark to light blue). The default value b = 2/3 is marked by 
the dotted line. 
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Figure 10. Equilibrium star formation efficiency per free fall time, eg, versus the turbulent production efficiencies ett (left panel), 
and egN (right panel) for the number densities n = 200, 100, 50, 37.5, 34.25, 25cm -3 (from dark to light blue). The default values 
(ett = 0.025, esN = 0.085) are marked by the dotted lines. 



18 H. Braun & W. Schmidt 



1 0.140 1- 
I 

^0.135 
I 

' 0.130h 



0.125 



50 100 150 

n [cm 3 ] 



200 



Figure 11. Relative deviation of the equilibrium solution for the 
warm-gas energy, M w ,eq, from the first two terms on the right- 
hand side of Eq. JS6P . corresponding to zero turbulent energy. 
The different curves are obtained for Ms = 0.0, 0.2, 0.4, 0.6, 0.8, 
1.0 from dark to light red. 
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Figure 13. Equilibrium star formation efficiency per free fall 
time, £ff , versus logarithmic specific energy of SN log(«sN/ u SN o)i 
where u S N,o = 6- 10 49 erg/M Q , for n = 25, 37.5, 50, 75, 100, 200 
cm -3 from light to dark blue. 
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Figure 12. Equilibrium star formation efficiency per free fall 
time, Eg, versus number density n for different values of the tur- 
bulence production efficiencies ett and esrj: ett = 0.025 (green), 
0.085 (black), 0.17 (red), and e S N = 0.025 (dotted), 0.085 (solid), 
0.17 (dashed). 

without external turbulence driving. This effect is caused by 
internal turbulence driving, which is considered next. 

Thermal instability and stellar feedback 

If the energy transfer from larger scales is small, internal 
driving by the thermal instability and stellar feedback dom- 
inate the production of turbulence. In this case, e t is mainly 
controlled by the turbulence production efficiencies ett and 
esN- An increase of ett lowers the thermal energy of the warm 
gas , u^, as the cooling instability transfers thermal energy 
to turbulence more efficiently. For higher esN, less energy is 
deposited by SNe in the warm phase and also more turbu- 
lent energy is produced, which tends to decrease the star 
formation rate. However, the effect of esN is limited for high 
ett because the increased production of turbulent energy by 
the thermal instability reduces star formation and, conse- 
quently, supernova feedback. On the other hand, only a rel- 
atively small fraction of the thermal energy of the warm 
phase can be converted into turbulent energy by the ther- 
mal instability without violating the second law of thermo- 
dynamics. This suggests that ett has to be small compared 
to unity. 

Since turbulence produced by SNe on length scales 
much smaller than I is rapidly dissipated into thermal en- 
ergy, only the fraction specified by esN effectively enters the 
turbulent energy e t , while the fraction 1 — esN is immediately 
turned into thermal energy. This implies that esN is scale- 
dependent. For I greater than a few parsec, esN must not be 
lower than a few percent. Otherwise the star formation effi- 



10 20 30 40 50 

Z[pc] 

Figure 14. Equilibrium star formation efficiency per free fall 
time, £ff , versus length scale £ for different number densities n = 
50, 75, 100, 200 cm~ 3 (from light to dark blue), without external 
driving. 

ciency would become significantly greater than 0.1, in con- 
tradiction to the m ajority of observations jKrumholz fc Tar] 
l2007l ; lMurravll201ll ). The average total energy deposited by 
a supernova in the interstellar medium is Ssn ~ 10 51 erg. 
Roughly 8.5 • 10 49 erg of this energy enter the ISM in form of 
kinetic energy, or fr om the perspective of our model in form 
of turbulent energy (|Thornton et al.ll 19981 ). Hence, we adopt 
the value esN = 0.085 as default. The left and middle panels 
of Fig. [10] show that eg saturates if one efficiency is much 
greater than the other. This is a consequence of p c , pa — > p, if 
turbulent energy is efficiently produced by whatever mech- 
anism. To obtain a plausible star formation efficiency, we 
choose e t t = 0.025. 

With the number of SNe per solar mass of formed stars, 
1 r^M Q dN 

n S N = — / dm,, (90) 

M * JsM e dm, 

we estimate the feedback energy per solar mass: 

u SN = n S NE SN ^-j^ « 6- 10 49 erg/M , (91) 

with the feedback fraction /3 as defined in Eq. [SS] This value 
is subject to large uncertainties. However, Fig. ll3l shows that 
the star formation efficiency is relatively robust against vari- 
ations in usn if its value is at least the same order of mag- 
nitude as the above default value for the model. 

Scale dependence 

For length scales I > 10 pc, the star formation efficiency 
is almost exactly scale- invariant (see Fig. 1141) . This demon- 
strates that even without maintaining a scale-invariant rate 
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Figure 15. Equilibrium star formation efficiency per free fall 
time, eg, versus the turbulence scaling exponent n for n = 
37.5, 50, 75, 100, 200 cm" 3 from dark to light blue. 



massive stars and the state of the absorbing gas. Both gas 
phases are affected by Lyc- heating (see Sect. 14. ip . and it 
reduces the amount of warm gas dropping into the cold 
phase if the cooling instability is active. As a consequence, 
an increase of £Lyc lowers the level of saturation of eg(n) 
significantly (see Fig. 1 1 T|) . Even for small x^yc, the den- 
sity threshold of star formation is affected by the influence 
of Lyc-heating on the thermal pressure of the warm phase 
and, thus, on p c , P a- Since the cross section decreases sig- 
nificantly for large excess energies above the Rydberg en- 
ergy and inelastic scattering processes become more likely, 
x~Lyc > 0.5 eV is not plausible. We use a;L yc = 0.1 eV as 
default value. 



of energy injection by external turbulence forcing, the sys- 
tem settles into an equilibrium state, for which the clump 
length scale l c and the Mach number of turbulence in the 
cold-gas phase, Ale, are nearly independent of I. Internal 
driving and turbulent dissipation regulate e t such that the 
scaling of the turbulent velocity fluctuations from I to l c by 
the power law (|12|l results in a fixed A4 C oc a c , regardless 
of the choice of I. However, the equilibrium is influenced by 
the scaling parameter rj. For weakly compressible turbulence 
(up to Mach numbers around unity), r\ is close to the Kol- 
mogorov value 1/3, whereas r\ rises to 1/2 for supersonic 
turbulence. The dependence of Eg on rj is shown in Fig. 1151 
Since we assume that turbulence is supersonic within the 
cold clumps, but transonic in the warm gas, and there are 
only little changes of Eg if rj is about 1/3, we set rj = 1/3 as 
a reasonable approximation if I is greater than l c . 

If I is only a few parsec or less, on the other hand, Z c may 
exceed I. In this case, basic assumptions of the model break 
down, as the notion of cold-gas clumps in pressure balance 
with the warm gas in the reference volume of size I becomes 
meaningless. The percolation of the cold phase and the tran- 
sition to a one-phase medium is not yet implemented in the 
model (see Sect. 15. ip . If the model in its present form is 
applied as an SGS model in AMR simulations, a maximum 
refinement limit has to be applied such that l c < I is ensured. 
In cosmological simulations with i m j n ~ 100 pc, this condi- 
tion will almost certainly be satisfied. For high-resolution 
simulations of individual galaxies with l m i n ~ 1 pc, however, 
situations, where cold-gas regions extend over several zones 
cannot be avoided. Then a viable model has to deal the 
transition to a one-phase medium. 

Photo-dissociation and heating 

The intensity of the ambient interstellar radiation field I v 
(in units of the Draine-field) is an external parameter and 
depends on the environment of the reference volume, i. e., 
the location and structure of the host galaxy and its sur- 
roundings. Thus, it is a parameter at the same level as n, Z 
and A4s, which has a direct impact on / c ,h 2 an d u w via HV 
dissociation and photoelectric heating. Fig. [16] demonstrate 
that the density threshold of star formation and the shape 
of Eg(n) change significantly with variations in I v , while 
the saturation level of eg remains nearly unaffected. Con- 
sequently, the threshold density of star formation is mainly 
determined by Z and I v . 

The heat deposited in the gas per Lyc-photon, x^ yG , de- 
pends on the shape of the spectrum emitted by the young 



Clump collisions 

The evaporation of clumps due to collisions lowers the p c /pw 
ratio, but the clumps are typically too small and collisions 
are too rare to influence the evolution significantly even for 
a high efficiency e cc . If l c becomes comparable to I, clump 
collisions cannot be applied for obvious reasons. We thus 
neglect this effect altogether by setting e cc = 0. 



Prestellar mass loss 

As mentioned by IChabrier fc Hennebelld l|201Ch . the CMF 
(the observed mass function of gravitationally bound cores 
in molecular clouds) and the IMF (the initial mass function 
of stars) are similar, except for an almost mass-independen t 
shift by factor about 2-3 (also see lMatzner fc McKeefeoOOh . 
We account for the mass reduction due to the evolution from 
the CMF to the IMF by reducing the star formation effi- 
ciency by a factor 1 — /i oss , where /i oss is interpreted as the 
mass fraction that is ejected during the collapse prior to star 
formation. As one can see in Fig. 1181 increasing /i oss reduces 
the star formation rate significantly. A good agreement with 
observational relations is obtained for the intermediate value 
/loss = 0.6. 



6.4 Comparison to observations 

Comparisons of one-zone results with observations are dif- 
ficult, mainly for the following reasons. Firstly, the con- 
version of the modelled volume densities into the corre- 
sponding surface densities is nontrivial. Secondly, star form- 
ing regions are generally not in local star formation equi- 
librium, which could explain the brea kdown of Kennicutt- 
Schmidt relations on small scales (e.g. lOnodera et alj|20ld : 
ISchruba et al ]|20ld : lMurravl201ll ). Without the detailed dy- 
namic environment in a hydrodynamic simulation, we can 
only draw conclusions on the basis of the equilibrium so- 
lutions calculated with our model. Even so, we are able to 
demonstrate that these solutions are consistent with the con- 
straints set by observations on kpc scales. 

Since the calculation of / c ,h 2 in °ur model only treats 
molecular hydrogen in cold clumps shielded from UV- 
radiation, we neglect the molecular hydrogen in radiation- 
dominated areas, where a certain amount exists in equilib- 
rium between production and radiative destruction. As a 
consequence, the predicted molecular hydrogen mass might 
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Figure 16. Equilibrium star formation efficiency per free fall time, eg , versus number density n for the normalised interstellar radiation 
field I v = 0.5, 1.0, 1.5, 2.0, 3.0 from dark to light blue (left panel). The metallicity of the gas is Z = O.^-Zq, and eg as function of log(/„) 
for n = 25, 50, 75, 100, 200 cm -3 from light to dark blue (right panel). The default I v = 1 is marked by the vertical dotted line. 
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Figure 17. Equilibrium star formation efficiency per free fall time, eg, versus number density n for x^ yc = 0.01, 0.1, 0.5, 1.0, 2.0, 4.0 
eV/photon from dark to light blue (left panel), and eg as function of log(x'Lyc/ e V) for different n = 25, 26.5, 28.125, 31.25, 37.5, 50, 
100 cm" 3 from light to dark blue (right panel). The default £Lyc = 0.1 eV is marked by the dotted line. 
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Figure 19. Molecular fraction /H 2 ,tot versus logarithmic total 
hydrogen surface number density log(27Vjj 2 + Nm) in star for- 
mation equilibrium for differe nt Z/Zp, = 1.0, 0.8, 0.6, .4 from 
dark to light blue. Red squ ares jRachford et al l2002l.l2009h, green 
circle s jWolfire et alj|2008h . and purple diamonds jGillmon et al.l 
I2006T) represent observations of translucent clouds in the Galaxy. 

be systematically too low. This discrepancy can become par- 
ticularly strong in the case of low (column) densities. To es- 
timate the hydrogen column density, we set (2Nh 2 + Nm) — 
Xpl/mn- Fig. 1191 shows that the transition from marginal 
to significant total molecular fractions /H 2 ,tot occurs at col- 
umn densities that are in good agreement with observations 
of translucent clouds in the Milky Way. 

Observed star formation surface densities and depletion 
timescales, as well as the related surface densities of molec- 
ular and atomic hydrogen, are projected quantities. In the 
simplest case of a face-on galactic disk, these quantities are 
averaged over the thickness of the disk and over a certain 
solid angle (or area). Since we can not account for spa- 
tial structures in the the equilibrium one-zone models, we 
plot of the star formation density over the atomic, molecu- 



lar, and total volume densities (with helium and metals in- 
cluded) for different metallicities and varying external driv- 
ing in Figs. 1201 1221 and !21l respectively. However, these plots 
should reproduce observed trends because the typical thick- 
ness of a star forming region is of the order of > 10 pc, which 
is comparable to the typical length scale £ of our models, and 
the gas in star forming regions contribute most to column 
densities. Ind eed, in comparison to observational rat es and 
densities (e.g. iBigiel et al.ll2008l ; ISchruba et al.ll201ll . Fig. 4 
and Fig. 11, respectively), we find a shift in numbers by a 
factor > 10 Q, in both directions, but the general behaviour 
is very similar. 

Figure [20l shows that the star formation rate is clearly 
not correlated to Hi-densities below < 1 Mqpc -3 (roughly 
corresponding to column densites < 10 Mqpc -2 ), and higher 
densities averaged over the cold and warm phases are atyp- 
ical. The star formation rate as a function of the total gas 
density (see Fig. I21[) switches from zero at low densities, for 
which the fraction of H2 is neglibile, to a tight correlation 
above ~ IMopc -3 . This threshold is caused by the tran- 
sition from atomic to shielded molecular gas. Correspond- 
ingly, the star formation rate mainly correlates with the H2- 
density, which agrees with the KMT09 model. For a given 
set of parameters, the molecular gas depletion timescale 
tde P ,H 2 = Ps/(/c,h 2 Pc) varies only little over two orders of 
magnitude in Eb-density (see Fig. l22p . While the model pre- 
dicts a depletion time of ~ 0.7 Gyr if internal driving is the 
dominant mechanism of turbulence production, the deple- 

4 The densitiy plots are in units of Mqpc -3 and M0Myr pc — 3 , 
while the observational column density plot are usually in units 
of M pc~ 2 and M0yr _1 kpc~ 3 . 
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Figure 18. Equilibrium star formation efficiency per free fall time, eg, versus number density n for /i oss = 0.0, 0.2, 0.4, 0.6, 0.8, 0.95 
from dark to light blue (left panel), and eg as function of /i oss for different n = 31.25, 37.5, 50, 75, 100cm -3 from light to dark blue 
(right panel). The default /i oaa = 0.5 is marked by the dotted line. 



tion time scale increas es significant l y for e xternal turbulence 
driving. For example, iBigiel et al.l (|201ll ) find ~ 2 Gyr ob- 
servationally, which could be maintained by external driving 
for H2-densities up to a few Mqpc -3 . However, as we un- 
derestimate the H2-content in the gas, the depletion time 
tends to be too low in our model toward lower densities so 
that less external energy injection might be required to ob- 
tain the observed depletion time. A further caveat is the 
assumption of equilibrium. 

Eventually, only the application as an SGS model in 
numerical simulations of galaxies will enable us to calculate 
relations between the star formation rate and the surface 
density by projecting the computed density fields. Since star 
forming regions go through different evolutionary stages, 
during which star formation occurs episodically, the star 
formation rate has to be integrated in time to incorporate 
non-equilibrium effects into the depletion time scale. 



7 DISCUSSION AND CONCLUSIONS 

In this paper we propose a model for the multi-phase 
ISM and star formation, considering the effects of tur- 
bulence and stellar f eedba ck. Based on the concept of 
ISpringel fc Hernquistl l|2005l ). we split the gas content of a 
region of given size into a two distinct fractions, represent- 
ing a diffuse warm and a clumpy cold component. However, 
our model goes significantly beyond their approach. By ap- 
plying a simplified treatment of molecular hydrogen forma- 
tion and destruction, we relate the star formation rate to 
the fra ctional density of mole cular hydrogen in the cold-gas 
phase |Krumholz et alj|2009h [KMT09] . While star forma- 
tion models that are applied in numerical simulations usu- 
ally assume a constant efficiency parameter that is globally 
calibrated against the Kennicutt- Schmidt law, we dynami- 
cally calculate the star formation effic iency on the basis of 
local p hysical processes. In the spirit of lKrumholz fc McKeei 
(2005), star formation is regulated by the virial parameter 
and the turbulent velocity dispersion of cold clumps. The in- 
terrelationship in our model is more elaborate though. Tur- 
bulent energy can be produced by a turbulent cascade from 
larger scale, but also via internal driving by the thermal 
instability of the gas and by supernova feedback. To deter- 
mine the gas fraction that collapses per free-fall timescale 
into stars, we assume a log-normal distribution of density 
fluctuations in the cold gas and we relate the critical den- 
sity for a gravitational collapse to the virial parameter and 



the turbulent M ach number on the typical length scale of 
the cold clumps (|Padoan fc Nordlun J201lh [PN11]. To ac- 
count for the effect of SNe, we use a delayed feedback model. 
Apart from the turbulent energy, the fractional densities of 
the cold and warm phases and the thermal energy of warm 
gas are evolved (the temperature of the cold gas is assumed 
to be constant). Mass and energy is exchanged between the 
phases via radiative cooling, heating and mixing processes. 
An important source of heating is the stellar population in 
the volume. We consider two feedback mechanisms, taking 
the time scales of stellar evolution into account: the Lyman- 
continuum emission and SNII-explosions of young massive 
stars. To close the system of equations, we assume an effec- 
tive (i. e., thermal plus turbulent) pressure balance between 
the cold clumps and the surrounding warm gas. 

By integrating the evolutionary equations for the aver- 
aged mass fractions and energies in a given spatial volume, 
we have obtained semi-analytical one-zone models, with the 
total gas density n, the metallicity Z and the rate of en- 
ergy injection by external turbulence forcing, £, as main 
parameters. Of particular interest are equilibrium solutions 
with a constant star formation rate. Fig. shows that, we 
obtain asymptotic Kennicutt-Schmidt-relations with slope 
1.5 (p B = pEff/tff oc p ) toward high densities, which is a 
consequence of the asymptotically constant star formation 
efficiency eg. Depending on the metallicity and other param- 
eters, the threshold densities are typically between 20 and 
about 200 cm -3 . For reasonable choices of the model coef- 
ficients that control internal turbulence driving and heat- 
ing, a star formation efficiency of around 0.5% is obtained 
above the t hreshold densities, in agreement with observed 
values (e.g. iKrumholz fc Tanll2007l; Bigiel et al.ll200Sl , l201ll ; 
ISchruba et al.ll201ll ; lOnodera et alj|20ich . External turbu- 
lence driving (i. e., energy transfer from larger scales via 
the turbulent cascade) decreases the star formation rate and 
slightly changes the slope of the power-law branches (Fig. 1211 
right panel). This is primarily caused by the effect of the 
turbulent pressure on the average density of the cold phase, 
while the direct influence of turbulence on the star forma- 
tion efficiency, following the prescription of PN11, plays a 
role in violently turbulent environments. In the latter case, 
also the production of molecular hydrogen fraction is af- 
fected via the turbulent clumping factor. Remarkably, the 
star formation efficiency is quite sensitive on the factor b in 
Eq. (|23l) for the width of the density pdf as a function of the 
tur bulent Mach number in the cold-gas phase. As shown 
by iFederrath et all J20ld ). b is related to the mixture of 
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Figure 20. Star formation rate as function of the fractional atomic density p — / Cj h 2 Pc for Z/Zq = 1.0, 0.8, 0.6, 0.4, 0.3 from dark to 
light blue (left panel) and for different turbulent energy injection rates E in terms of My = 0.0, 0.4, 0.8, 1.2, 1.4 (from dark to light 
blue, right panel). The long- and short-dashed lines indicate the asymptotes with slopes 1.4 (My = 0.0, all Z/Zq) and 1.6 {My: = 1-0), 
respectively. The orange dot-dashed lines mark fixed depletion time scales of 0.1, 1, and 10 Gyr from the top to the bottom of the graph. 
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Figure 21. Star formation rate as function of the total density p for Z/Zq = 1.0, 0.8, 0.6, 0.4, 0.3 from dark to light blue (left panel) 
and for different turbulent energy injection rates £ in terms of My = 0.0, 0.4, 0.8, 1.2, 1.4 (from dark to light blue, right panel). The 
Kennicutt-Schmidt-asymptote with slope 1.5 is indicated by the dashed line, orange dot-dashed lines as in Fig. 1201 



the solenoidal and compressive components of the turbulent 
velocity field. Moreover, they concluded from comparisons 
with observed two-point statistics of turbulence in molecu- 
lar clouds that this mixture varies for different clouds. Thus, 
it appears to be important to account for variations in the 
turbulence statistics. 

Recent observations indicate a particularly tight cor- 
relation of the star formation rate with the molecular gas 
column densities in galaxies down to kpc scales (KMT09 
give an overview of observational results) . Since we consider 
local regions of the ISM with a size smaller than the galactic 
disc thickness, it is not reasonable to express the results from 
our one-zone models in terms of column de nsities. For the 
same reason, comparisons with the model of lOstriker et al.1 
(2010) are difficult. Nevertheless, we find that the relation 
between the equilibrium star formation rate and the density 
of molecular hydrogen closely follows a power law, particu- 
larly for solar metallicity. As one can see in the left panel of 
Fig. [22}, the star formation rate is p B tx p 1 for sufficiently 
high density. Strong external turbulence forcing significantly 
reduces the star formation rate and the slope of the asymp- 
tote increases from 1.4 to about 1.6 (Fig. [22] right panel). 
In this regard, it is interesting that KMT09 distinguish two 
different regimes, in which molecular clouds are either self- 
regulated (at low surface densities) or significantly affected 
by their galactic environment (at high surface densities). In 
the former case, they derive a slope of about 1.4, whereas the 
slope is about 1.6 in the latter case. A plausible interpreta- 
tion in the context of our model is that these regimes roughly 
correspond to internal turbulence driving as the dominating 
production mechanism (negligible E) vs. significant turbu- 



lence production by the transport from instabilities on large 
scale to molecular cloud scales (large E). To corroborate 
this interpretation, the model has to be applied in simula- 
tions of disc galaxies. The modelled equilibrium star forma- 
tion rates and depletion time scales a re roughly consisten t 
with those found observationally (e.g. ISchruba et al1l201ll ). 
and the modelled relation between star formation rate and 
molecular gas density is in agreement with a more or less 
constant molecula r gas depletion time scale as observed by 
iBigiel et al.l (|201ll ). 

In such simulations as well as in cosmological simula- 
tions, it is common to assume a constant star formation effi- 
ciency beyond a certain threshold density. The results of our 
numerical study suggest that this is a reasonable approxi- 
mation. However, rather than using this as an entirely phe- 
nomenological input to the simulations, the equilibrium val- 
ues of the star formation efficiency calculated with our model 
follow from the sub-resolution physics of the ISM. Moreover, 
rather than setting stiff density thresholds, the model yields 
transition values depending on the varying gas density in nu- 
merical simulations. To utilise the equilibrium solutions as 
a parametrization of star formation, tables of the star for- 
mation efficiency as function of density and metallicity can 
be calculated (the code calculating the efficiencies can be an 
be obtained from the authors upon request). The rate of ex- 
ternal turbulent energy production could be estimated, for 
inst ance, from the large- scale velocity dispersion in galaxies 
(see iBurkert et all |2010| ) . 

While such a simplified approach has its merits, it can- 
not account for dynamical effects. A crucial problem is the 
calculation of the local rate of turbulent energy production 
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Figure 22. Star formation rate as function of the fractional molecular density f c ,n 2 Pc for Z/Zq = 1.0, 0.8, 0.6, 0.4, 0.3 from dark to 
light blue (left panel) and for different turbulent energy injection rates E in terms of A^s = 0.0, 0.4, 0.8, 1.2, 1.4 (from dark to light 
blue, right panel). The long- and short-dashed lines indicate the asymptotes with slopes 1.4 = 0.0, all Z/Zq) and 1.6 (My: = 1-0), 

respectively. Orange dot-dashed lines are as in Fig. 1201 



on the grid scale due the shear of numerically resolved tur- 
bulent flow in a simulation (i. e., the energy transfer from 
length scale greater than the size of the grid cells to unre- 
solved length scales). This is the meaning of E if the pro- 
posed model is applied as a sub-grid scale model. Since the 
turbulent velocity fluctuations in the ISM can assume a sig- 
nificant fraction of the sound speed or even become super- 
sonic, an incompressible tur bulence model is not sufficient. 
ISchmidt fc Federrathl (|201ll ) [SF11] provide a formula for E 
in the highly compressible regime. A complete model for the 
turbulent multi-phase ISM and star formation is obtained by 
rewriting Eqs. (J42] [61] [62j [73] [75] [46j) as partial differen- 
tial equations with fluid-dynamical advection terms, where 
the length scale I is given by size of the grid cells, A, and 
e t = e sgs is identified with the unresolved fraction of the 
kinetic energy (see SF11). These equations supplement the 
Euler equations for the total gas density, the momentum, 
and the total energy. Solving the complete set of equations 
will be a substantial numerical challenge. 

To perform simulations of galaxies in cosmological en- 
viron ments, adap t ive m esh refinement (AMR) is indispens- 
able. |m 

aier et all (|2009l ) incorporated an SGS turbulent en- 
ergy equation for moderately compressible turbulence into 
AMR simulations of galaxy clusters. This method also can 
be applied using the multi-phase model for the turbulent 
ISM. Then the length scale I of the model corresponds to 
the varying grid scale, and the scale-dependent turbulent 
energy has to be adjusted if refined grids are inserted or 
solutions on finer grids are projected to coarser grid levels. 

An advantage of our model is that components can be 
modified, replaced and added as our understanding of the 
physics of the ISM progresses. For example, the computa- 
tion of the dimensionless star formation rate from the tur- 
bulent cold-gas density pdf (see Sect. 13. 1 ^ is more or less 
heuristic. We anticipate that theoretical advances and re- 
sults form small-scale simulations will soon lead to improve- 
ments. An important issue we have not considered so far is 
the influence of magnetic fields. The role of MHD turbulence 
is already emphasized by PN11. Although additional com- 
plications arise when magneto-turbulent fluctuations have 
to be treated on sub-grid scales, it is a problem that can be 
tackled. Furthermore, a more d etailed treatment of chemica l 
processes is desirable, although iKrumholz fc Gnedinl (|201ll ) 
have demonstrated that the simple analytical model for the 
molecular hydrogen fraction in KMT09 agrees quite well 
with an explicit reaction network in cosmological simula- 



tions, at least if the metallicity is not much lower than solar. 
The multi-phase model also has to be adapted to the sim- 
ulation framework. For cosmological simulations with rela- 
tively coarse resolutions, the neutral gas phases should be 
embedded in a hot ionised medium of low density. The other 
extreme are simulations of isolated disc galaxies with very 
high resolution, in which cold clumps can be marginally re- 
solved so that several neighbouring grid cells are completely 
filled by cold gas. In this case, it is necessary to switch from 
the two-phase description to the one- phase limit. 

The predictive power of astrophysical simulations, in 
which the ISM is only partially resolved, will increase by 
applying the equilibrium solutions or by implementing the 
full multi-phase SGS model. This will allow us, in turn, to 
test and to modify the underlying physical assumptions. 
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